Loading libraries

Loading Data

Reviewed May 15th, 2025:

Plots for Gross Photosynthesis - Cabral:

  • Cabral is significant for silicate at linear (0.04) and nonlinear (0.02)
  • Cabral is only significant for SGD Dils in the nonlinear model (0.02)
  • With removing an outlier point in the last dilution, Cabral is significant for GP both linear (0.015) and nonlinear (0.015)
# mutate(ammonia_umolL = ifelse(ammonia_umolL = "<0.02", 0.01, ammonia_umolL)) %>% 

###################
## Silicate
####################

### confirm linear model 
GP_Cabral_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  #filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(silicate_umolL + 1) 0.060043 0.060043     1 110.84  3.2591 0.07374 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm nonlinear model 
GP_Cabral_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
 # filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.20901 0.10451     2 110.86  6.1685 0.002883
##                                    
## poly(log(silicate_umolL + 1), 2) **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### both are significant so can add either trend line 

GP_Cabral_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 # geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log silicate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with Silicate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_Silicate

###################
## SGD dils
####################

### confirm a linear model 
GP_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.045966 0.045966     1 46.338   3.266 0.07722 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
GP_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)  
## poly(log(SGD_number + 1), 2) 0.10382 0.051908     2 45.272  3.9695 0.0258 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GP_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_SGDDils

###################
## SGD dils WITHOUT POTENTIAL OUTLIERS
####################
GP_Cabral_SGDDils_nooutlier <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
 scale_color_manual(values=pnw_palette("Bay", n=9)) +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  #geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
  geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_SGDDils_nooutlier

### confirm a linear model 
GP_Cabral_SGD_linearmodel_nooutlier <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_linearmodel_nooutlier)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.081847 0.081847     1    50  6.2821 0.01549 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
GP_Cabral_SGD_nonlinearmodel_nooutlier <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_nonlinearmodel_nooutlier)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## poly(log(SGD_number + 1), 2) 0.11464 0.057321     2    49  4.5402 0.01553 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plots for Gross Photosynthesis - Varari

  • Varari is ONLY significant for silicate in nonlinear model (0.02)
  • Varari has no significant relationship with SGD and GP
###################
## Silicate
####################

### confirm linear model 
GP_Varari_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>%
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## log(silicate_umolL + 1) 0.015463 0.015463     1 62.55  0.7003 0.4059
### confirm nonlinear model 
GP_Varari_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.15527 0.077633     2 61.577   3.925 0.02487
##                                   
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### ONLY significant in nonlinear model 

GP_Varari_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log silicate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Varari with Silicate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Varari_Silicate

###################
## SGD dils
####################

### confirm a linear model 
GP_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.010297 0.010297     1 60.065  0.4631 0.4988
### confirm a nonlinear model 
GP_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.089681 0.044841     2 59.102  2.1119   0.13
GP_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
#  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Varari_SGDDils

Plots for Respiration - Cabral:

  • There is no significance at Cabral for any environmental parameter or broad SGD dilutions
###################
## SGD dils
####################

### confirm a linear model 
R_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 9.037e-08 9.037e-08     1 46.085       0 0.9962
### confirm a nonlinear model 
R_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.011554 0.0057771     2 45.076  1.5607 0.2211
R_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "R Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
R_Cabral_SGDDils

Plots for Respiration - Varari

  • There is no significance at Varari for any environmental parameter or broad SGD dilutions
###################
## SGD dils
####################

### confirm a linear model 
R_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.0066251 0.0066251     1 59.962  2.1962 0.1436
### confirm a nonlinear model 
R_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.012828 0.0064138     2 58.984  2.1691 0.1233
R_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
#  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "R Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
R_Varari_SGDDils

Plots for Calcification - Cabral:

  • Cabral was allegedly significant for temperature (0.01) and pH (0.01) with a nonlinear relationship - but not seeing that here for temp, only for pH
  • Cabral is NOT significant for SGD Dils
###################
## pH
####################

### confirm nonlinear model 
C_Cabral_pH_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(new_pH+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_pH_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## poly(log(new_pH + 1), 2) 0.043887 0.021943     2 47.288  3.2971 0.04566 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Cabral_pH <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=new_pH, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log pH",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Cabral with pH") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Cabral_pH

###################
## SGD dils
####################

### confirm a linear model 
C_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.015433 0.015433     1 46.137  2.1383 0.1504
### confirm a nonlinear model 
C_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.01618 0.0080902     2 45.126  1.0993 0.3418
C_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  #geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Cabral_SGDDils

Plots for Calcification - Varari

  • Varari has a significant linear and nonlinear relationship with phosphate (<0.0001) and (0.0002) and TA (0.001)
  • Varari has a significant … relationship with SGD
###################
## Phosphate
####################

### confirm nonlinear model 
C_Varari_phosphate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(phosphate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
    filter(phosphate_umolL<1) %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_phosphate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                     Sum Sq   Mean Sq NumDF  DenDF F value
## poly(log(phosphate_umolL + 1), 2) 0.011868 0.0059339     2 58.832  1.6788
##                                   Pr(>F)
## poly(log(phosphate_umolL + 1), 2) 0.1954
C_Varari_Phosphate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
   filter(phosphate_umolL<1) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=phosphate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log phosphate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Varari with Phosphate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_Phosphate

###################
## TA 
####################

### confirm a nonlinear model 
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(TA_initial+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
    filter(TA_initial<2550) %>% #### IF CHOOSE TO REMOVE HIGH POINTS, ONLY SIGNIFICANT IN NONLINEAR (otherwise significant in linear too)
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq  Mean Sq NumDF  DenDF F value   Pr(>F)   
## poly(log(TA_initial + 1), 2) 0.056876 0.028438     2 59.115  5.1362 0.008776 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_TA <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  # filter(TA_initial<2550) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=TA_initial+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "log TA",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Varari with TA (Minus the potential outliers)") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_TA

###################
## Silicate 
####################

### nonlinear model 
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.057443 0.028721     2 60.29  4.7675 0.01196
##                                   
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## not poly 
C_Varari_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)
## log(silicate_umolL + 1)      0.046571 0.046571     1 107.262  6.9107 0.009826
## site                         0.000212 0.000212     1  47.521  0.0315 0.859918
## log(silicate_umolL + 1):site 0.002363 0.002363     1 107.262  0.3507 0.554987
##                                
## log(silicate_umolL + 1)      **
## site                           
## log(silicate_umolL + 1):site   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  # filter(TA_initial<2550) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=site)) +
  scale_color_manual(values=pnw_palette("Bay", n=2), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "log Silicate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Varari with Silicate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_silicate

###################
## SGD dils
####################

### confirm a linear model 
C_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.041992 0.041992     1 60.038  6.7823 0.01159 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>%
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq  Mean Sq NumDF  DenDF F value   Pr(>F)   
## poly(log(SGD_number + 1), 2) 0.082704 0.041352     2 59.049  7.3684 0.001391 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "NEC Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_SGDDils

Plot nutrient levels with SGD dils to see how different the nutrients actually are across the dilutions

  • just using as a check/baseline

Plot how many corals are net calcifying vs net dissolving

  • this is supplementarl figure used as additional evidence to show difference in Varari gw vs Cabral
  • Results: more fragments in Varari were net dissolving compared to Cabral (and at lower dilutions)
### should make a stacked bar plot with all corals in one group of dilution - % of which of those are positive or negative 

###################################
## attempt 2 
###################################
RespoR_Normalized_Full_withNEC_nutrients_calcifplot <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>%
  select(!c("phosphate_umolL":"TA_initial", "salinity", "date")) %>%
  mutate(rate_category = ifelse(umol.cm2.hr > 0, "positive", "negative")) %>% 
  group_by(SGD_number, site) %>% 
  dplyr::summarize(positive_count = sum(rate_category == "positive"),
            negative_count = sum(rate_category == "negative"))

RespoR_Normalized_Full_withNEC_nutrients_calcifplot2 <- RespoR_Normalized_Full_withNEC_nutrients_calcifplot %>% 
  pivot_longer(cols = c(positive_count, negative_count),
               names_to = "rate_category",
               values_to = "count") %>% 
  mutate(total_count = sum(count),
         proportion = count / total_count)

## plot 
Net_CalcvsDissolution2 <- RespoR_Normalized_Full_withNEC_nutrients_calcifplot2 %>% 
  ggplot(aes(x=as.factor(SGD_number), 
             y=proportion,
             fill=rate_category)) +
  geom_bar(stat = "identity", position = "fill", width=0.4) +
  #scale_fill_manual(values =  ifelse (umol.cm2.hr > 0, "red", "black")) +
  scale_fill_manual(values=c("red", "black"), 
                    labels = c("positive_count" = "Net Calcifying", "negative_count" = "Net Dissolving")) +
  theme_bw() + 
 # coord_trans(x="log") +
 # geom_hline(yintercept = 0) +
  facet_wrap(~site) +
 #coord_trans(x ="log") +
  labs(x = "% SGD by Volume",
       y = "Proportion of Corals",
      # title = "Net Calcifying vs Net Dissolving Corals per SGD Dilution", 
       fill = "Calcification State") +
   theme(axis.text.x=element_text(size=14), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20)) 
Net_CalcvsDissolution2

ggsave(here::here("Outputs", "Results", "Final", "NetCalcifying_NetDissolving.jpg"), 
       width = 10, height=7)

########################
## statistics 
#######################
CalcVsDissol <- lm(proportion ~ SGD_number*site, data = RespoR_Normalized_Full_withNEC_nutrients_calcifplot2)

anova(CalcVsDissol)
## Analysis of Variance Table
## 
## Response: proportion
##                 Df  Sum Sq  Mean Sq F value Pr(>F)
## SGD_number       1 0.00000 0.000000  0.0000 1.0000
## site             1 0.03840 0.038405  1.0651 0.3098
## SGD_number:site  1 0.00068 0.000684  0.0190 0.8913
## Residuals       32 1.15385 0.036058
########################
## calculate percentage that are dissolving 
#######################



percentDissolve <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         site=="Cabral") %>% 
  dplyr::mutate(percent_dissolved = (sum(umol.cm2.hr < 0) / n()) * 100)
  




##############################################
## look at raw data to see why dissolving is occurring at 0.01 
##############################################
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="0.01") ### looks like it is only 3 colonies, 1 at Cabral (Col 5, but barely at -0.02), and 2 at Varari (Col 3 and 8)

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="0.03") ### only the same 2 colonies at Varari (Col 3 and 8)

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="0.05") ### only Varari Col 3 

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="0.1") ### only Varari Col 3 

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="0.5") ### Varari Col 3 and barely Col 8 (0.02)

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="1")  ### Varari Col 1 and barely Col 3

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="2")  ### one colony at Cabral (Col 10), Varari Col 3 and Col 8

RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C", 
         SGD_number=="4")  ### one colony at Cabral (Col 8), Varari Cols 2,3 and 4

###################################
## check GP rates of Varari Colony 3  
###################################
RespoR_Normalized_Full_withNEC_nutrients_GPrawdata <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="GP", 
         new_colonynumber=="3"| new_colonynumber=="2")

RespoR_Normalized_Full_withNEC_nutrientsVGPonly <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="GP", 
         site=="Varari") %>% 
   filter(!is.infinite(umol.cm2.hr))

max(RespoR_Normalized_Full_withNEC_nutrientsVGPonly$umol.cm2.hr)
## [1] 1.400944
min(RespoR_Normalized_Full_withNEC_nutrientsVGPonly$umol.cm2.hr)
## [1] 0.1810457
###################################
## attempt 1 
###################################

Net_CalcvsDissolution <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  mutate(SGD_number = as.character(SGD_number), 
         rate_category = ifelse(umol.cm2.hr > 0, "positive", "negative")) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number, 
             y=umol.cm2.hr,
             fill=rate_category)) + 
  geom_col(width = 0.8,  colour = NA) + 
  #scale_fill_manual(values =  ifelse (umol.cm2.hr > 0, "red", "black")) +
  scale_fill_manual(values=c(positive='black',negative='red')) +
  theme_bw() + 
  geom_hline(yintercept = 0) +
  facet_wrap(~site) +
 #coord_trans(x ="log") +
  labs(x = "SGD values",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "Net Calcifying vs Net Dissolving Corals per SGD Dilution") +
   theme(axis.text.x=element_text(size=14), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
Net_CalcvsDissolution

Plot GP:R as a function of SGD

  • Cabral is significant (0.025) but Varari is not
## First make new df
GP_RRatio_bySGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  dplyr::select(!"sample_ID") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  group_by(P_R, SGD_number, date, site, new_colonynumber) %>% 
  filter(P_R!="C", P_R!="NP") %>% 
  #summarize(GP_R = GP/R)
  pivot_wider(names_from = P_R, values_from = umol.cm2.hr) %>% 
  mutate(GP_R = GP/R) 

### plots 
GP_R_bySGD_bothsites <- GP_RRatio_bySGDDils %>% 
  #filter(site=="Varari") %>% 
  mutate(new_colonynumber= as.factor(new_colonynumber)) %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=GP_R)) + 
  geom_point(size=3, aes(color=new_colonynumber)) +
  theme_bw() + 
  coord_trans(x="log") +
  facet_wrap(~site, scales = "free_y") +
  geom_smooth(method="lm", formula = "y~log(x)") +
  labs(x = "SGD Dilutions",
       y = "Gross Photosynthesis/Respiration Rates",
       title = "SGD by GP:R Relationship for Varari") +
   theme(axis.text.x=element_text(size=12), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_R_bySGD_bothsites

## model testing 
model_GP_R_SGD_Varari <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
                                        filter(site=="Varari"))
anova(model_GP_R_SGD_Varari)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.26869 0.26869     1 60.003  3.6316 0.06148 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_GP_R_SGD_Cabral <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
                                        filter(site=="Cabral"))
anova(model_GP_R_SGD_Cabral)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
## log(SGD_number + 1) 0.28731 0.28731     1 46.031  5.3534 0.0252 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NEW UPDATES FROM JUNE 19TH (MEETING WITH NYSSA):

Three main figures for publication will include:

Figure 1: Conceptual (made in Biorender)

Figure 2: ReVamp’d (correlation matrix) (ignore these)
  • this has all the significant correlations and labeled for each site
###################################
## Varari 
###################################

## data 
#corr_matrix_Varari <- RespoR_Normalized_Full_withNEC_nutrients %>% 
 # mutate(N_P = NN_umolL/phosphate_umolL) %>% 
 # filter(P_R=="R") %>% 
 # mutate("Salinity"=salinity, 
      #   "TA" = TA_initial, 
      #   "pH" = new_pH, 
      #   "N+N" = NN_umolL, 
       #  "Phosphate" = phosphate_umolL,
       #  "Silicate" = silicate_umolL) %>% 
         #"N+N:P" = N_P) %>% 
 # filter(site=="Varari") %>% 
 # select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", #"umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) 

# Compute a correlation matrix
#corr_Varari <- round(cor(corr_matrix_Varari, 
                  #       method="spearman"), 1)

# Compute a matrix of correlation p-values
#p.mat_Varari <- cor_pmat(corr_Varari)

# Visualize the correlation matrix
#ggcorrplot(corr_Varari, method = "circle")

# Reordering the correlation matrix using hierarchical clustering
#ggcorrplot(corr_Varari, hc.order = TRUE, outline.col = "white", method = "circle")

# Types of correlogram layout - Get the lower triangle
#ggcorrplot(corr_Varari, hc.order = TRUE, type = "lower", outline.col = "white", method = "circle")

# Change colors and theme
# Argument colors
#Varari_corr_matrix_plot <- ggcorrplot(corr_Varari, 
                                #      hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
                                  #    type = "lower", ## correlogram layout - for the lower triangle
                                   #   outline.col = "white",
                                   #   # method = "circle", ## determines the shape
                                  #    ggtheme = ggplot2::theme_classic, ## changes the theme 
                                  #    colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
                                   #   sig.level = 0.05, ## default p-value
                                   #   lab = TRUE, ## adds correlation coefficient
                                   #   p.mat = p.mat_Varari, ## adds significance
                                    #  insig = "blank", ## makes insig values blank
                                    #  title = "Pearson Correlation Matrix for Environmental Parameters at Varari",
                                #      legend.title = "Correlation", 
                                 #     tl.srt = 45, ## slant of x axis 
                                 #     tl.cex = 15, ## size of text axes  
                                 #     lab_size = 5)  + 
  #theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
      #  legend.title = element_text(size = 10, face="bold"),
      #  axis.text.x = element_text(size = 12, face="bold", color="black"),
      #  axis.text.y = element_text(size = 12, face="bold", color="black"))

#Varari_corr_matrix_plot
#ggsave(here::here("Outputs", "Results", "REDO_VarariCorrelation.jpg"), 
     #  width=10, height=7)


###################################
## Cabral 
###################################

## data 
#corr_matrix_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>% 
 # mutate(N_P = NN_umolL/phosphate_umolL) %>% 
 # filter(P_R=="R") %>% 
#  mutate("Salinity"=salinity, 
      #   "TA" = TA_initial, 
      #   "pH" = new_pH, 
       #  "N+N" = NN_umolL, 
        #  "Phosphate" = phosphate_umolL,
       #   "Silicate" = silicate_umolL) %>%  
        # "N:P" = N_P) %>% 
 #  filter(site=="Cabral") %>% 
 #  select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", # "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) 

# Compute a correlation matrix
# corr_Cabral <- round(cor(corr_matrix_Cabral, 
                  #        method="spearman"), 1)


# Compute a matrix of correlation p-values
# p.mat_Cabral <- cor_pmat(corr_Cabral)

# Visualize the correlation matrix
# ggcorrplot(corr_Cabral, method = "circle")
# 
# Change colors and theme
# Argument colors
# Cabral_corr_matrix_plot <- ggcorrplot(corr_Cabral, 
                                 #      hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
                                #       type = "lower", ## correlogram layout - for the lower triangle
                                 #      outline.col = "white",
                                      # method = "circle", ## determines the shape
                                 #      ggtheme = ggplot2::theme_classic, ## changes the theme 
                                 #      colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
                                      #sig.level = 0.05, ## default p-value
                                  #     lab = TRUE, ## adds correlation coefficient
                                  #     p.mat = p.mat_Cabral, ## adds significance
                                    #  insig = "blank", ## makes insig values blank
                                  #     title = "Pearson Correlation Matrix for Environmental Parameters at Cabral",
                                  #     legend.title = "Correlation", 
                                  #     tl.srt = 45, ## slant of x axis 
                                  #     tl.cex = 15, ## size of text axes  
                                   #    lab_size = 5)  + 
 #  theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
 # #        legend.title = element_text(size = 10, face="bold"),
      #   axis.text.x = element_text(size = 12, face="bold", color="black"),
      #   axis.text.y = element_text(size = 12, face="bold", color="black")) + 
 #  labs(x="Location", 
   #     y = expression(Phosphate~(mu*mol~L^-1))) 

# Cabral_corr_matrix_plot
# ggsave(here::here("Outputs", "Results", "REDO_CabralCorrelation.jpg"), 
       # width=10, height=7)
Figure 2: ReVamp’d again (ignore these)
  • this is revamped but still includes all the squares further interpretation
 ###################################
## Varari 
###################################

## data from above 

# Compute a correlation matrix
# corr_Varari2 <- round(cor(corr_matrix_Varari), 1)

# Compute a matrix of correlation p-values
# p.mat_Varari2 <- cor_pmat(corr_Varari2)

# Specify the desired order of parameters
#param_orderV <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names

# Reorder the correlation matrix
#corr_Varari <- corr_Varari[param_orderV, param_orderV]

# Visualize the correlation matrix
# ggcorrplot(corr_Varari, method = "circle")


# Varari_corr_matrix_plot2 <- ggcorrplot(corr_Varari, 
                                     # hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
                                  #     type = "lower", ## correlogram layout - for the lower triangle
                                  #     outline.col = "white",
                                      # method = "circle", ## determines the shape
                                  #     ggtheme = ggplot2::theme_classic, ## changes the theme 
                                  #     colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
                                      #sig.level = 0.05, ## default p-value
                                  #     lab = TRUE, ## adds correlation coefficient
                                     #p.mat = p.mat_Varari, ## puts an X on the insignificant ones 
                                    #  insig = "pch", ## makes insig values blank
                                     #insig = label_sig, 
                                  #    title = "Pearson Correlation Matrix for Environmental Parameters at Varari",
                                  #     legend.title = "Correlation", 
                                  #     tl.srt = 45, ## slant of x axis 
                                  # # #     tl.cex = 15, ## size of text axes  
                                   #    lab_size = 5) +
  #labs(y="Varari") +
 #  theme(legend.title = element_text(size = 10, face="bold"),
   #      axis.text.x = element_text(size = 12, face="bold", color="black"),
   #      axis.text.y = element_text(size = 12, face="bold", color="black"))

# Varari_corr_matrix_plot2
# ggsave(here::here("Outputs", "Results", "REDO_VarariCorrelation_extended.jpg"), 
     #   width=10, height=7)
# # 

###################################
## Cabral 
###################################

# Compute a correlation matrix
# corr_Cabral2 <- round(cor(corr_matrix_Cabral), 1)

# Compute a matrix of correlation p-values
# p.mat_Cabral2 <- cor_pmat(corr_Cabral2)

# convert to long format so can add asterisk 
#cor_long_Cabral <- melt(corr_Cabral2)
#p_long_Cabral <- melt(p.mat_Cabral2)

# Mark significant correlations with asterisks
#cor_long_Cabral$signif <- ifelse(p_long_Cabral$value < 0.05, "bold", "plain")

# Specify the desired order of parameters
# param_orderC <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names

# Reorder the correlation matrix
# corr_Cabral2 <- corr_Cabral2[param_orderC, param_orderC]

# ggcorrplot(corr_Cabral2, method = "circle")


# Change colors and theme
# Argument colors
# Cabral_corr_matrix_plot2 <- ggcorrplot(corr_Cabral2, 
                                     # hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
                                     #  type = "lower", ## correlogram layout - for the lower triangle
                                     #  outline.col = "white",
                                      # method = "circle", ## determines the shape
                                    #   ggtheme = ggplot2::theme_classic, ## changes the theme 
                                  #     colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
                                      #sig.level = 0.05, ## default p-value
                             #         lab = TRUE, ## adds correlation coefficient
                                     # p.mat = p.mat_Cabral2, ## adds significance
                                      #insig = "blank", ## makes insig values blank
                                     # title = "Pearson Correlation Matrix for Environmental Parameters at Cabral",
                                   #    legend.title = "Correlation", 
                                   #    tl.srt = 45, ## slant of x axis 
                                    #   tl.cex = 15, ## size of text axes  
                                    #   lab_size = 5)  + 
  # theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"), 
       #  legend.title = element_text(size = 10, face="bold"),
       #  axis.text.x = element_text(size = 12, face="bold", color="black"),
       #  axis.text.y = element_text(size = 12, face="bold", color="black")) #+ 
  #geom_text(data = cor_long_Cabral, aes(x = Var1, y = Var2, label = signif), color = "black", size = 5, inherit.aes = FALSE)


# Cabral_corr_matrix_plot2
# ggsave(here::here("Outputs", "Results", "REDO_CabralCorrelation_extended.jpg"), 
      #  width=10, height=7)


########################### 
## patch together 
###########################

#### remove x-axis labels 
# p_Ccorr <- Cabral_corr_matrix_plot2 
# p_Vcorr <- Varari_corr_matrix_plot2 
 
# combined_corr_plot <- p_Ccorr/p_Vcorr + plot_layout(guides = 'collect')

# combined_corr_plot + 
#   plot_annotation(tag_levels = 'A') &
#   theme(plot.margin = unit(c(1, 1, 1, 1), "cm"), 
 #        plot.tag = element_text(size = 18, face = "bold"))

# ggsave(here::here("Outputs", "Results", "Final", "CorrMatrices.jpg"), 
     #   width=12, height=15)

Figure 2a: correlation matrix ONLY with silicate AND spearman correlation - THIS IS THE ONE TO USE

  • this is the final figure included for publication, which only has silicate compared to other env params AND uses the Spearman correlation which is better for nonlinear, monotonic data
 ###################################
## Varari 
###################################
## data 
corr_matrix_Varari_silicateonly <- RespoR_Normalized_Full_withNEC_nutrients %>% 
#  mutate(salinity = if_else(date == "6/17/23", salinity - 1, salinity), ### correcting for the weird day of salinity on 6/17 that was way higher than any other days. See code in Salinity_pH_TA_TestingRelationships for the deduction/evidence. Edited on 03/21/2025
  mutate(N_P = NN_umolL/phosphate_umolL) %>% 
  filter(P_R=="C") %>% 
  mutate("Salinity"=salinity, 
         "TA" = TA_initial, 
         "pH" = new_pH, 
         "N+N" = NN_umolL, 
         "Phosphate" = phosphate_umolL,
         "Silicate" = silicate_umolL) %>% 
         #"N+N:P" = N_P) %>% 
  filter(site=="Varari") %>% 
  select(!c("salinity", "silicate_umolL", "log_silicate", "log_silicate2", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) 

## version 2.0 
corr_Varari_spearman <- round(cor(corr_matrix_Varari_silicateonly, method="spearman"), 1)

# Compute a matrix of correlation p-values
p.mat_Varari_spearman <- cor_pmat(corr_Varari_spearman)

# Extract the correlations involving "Silicate"
corr_silicate_V <- corr_Varari_spearman[, "Silicate", drop = FALSE]

# Transpose the matrix to make "Silicate" the row instead of the column
corr_silicate_matrix_V <- t(as.matrix(corr_silicate_V))

# Convert the matrix to a long format
corr_melted_V <- melt(corr_silicate_matrix_V)

# Plot the correlation matrix vertically
Varari_corr_matrix_spearman <- ggplot(corr_melted_V, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile(color = "white") +  # Adjust width of the tiles to make the plot narrower
  scale_fill_gradient2(low = "steelblue3", high = "tomato3", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Spearman\nCorrelation") +
  geom_text(aes(label = round(value, 2)), color = "black", size = 5) +
  theme_minimal() + 
  theme(axis.text.x = element_text(vjust = 0.5, size = 12, face = "bold", color="black"),
        axis.text.y = element_text(size = 12, face = "bold", color="black")) +  
  coord_flip() +  # Flip coordinates for vertical layout
  labs(x = "", y = "")

# Display the plot
Varari_corr_matrix_spearman

ggsave(here::here("Outputs", "Results", "Final", "Vsilicate_corrmatrix.jpg"), 
       width=7, height=12)

###################################
## Cabral 
###################################
## data 
corr_matrix_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  mutate(N_P = NN_umolL/phosphate_umolL) %>% 
  filter(P_R=="C") %>% 
  mutate("Salinity"=salinity, 
         "TA" = TA_initial, 
         "pH" = new_pH, 
         "N+N" = NN_umolL, 
         "Phosphate" = phosphate_umolL,
         "Silicate" = silicate_umolL) %>%  
        # "N:P" = N_P) %>% 
  filter(site=="Cabral") %>% 
  select(!c("salinity", "log_silicate", "log_silicate2", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) 

# Compute a correlation matrix
corr_Cabral_spearman <- round(cor(corr_matrix_Cabral, method="spearman"), 1)

# Compute a matrix of correlation p-values
p.mat_Cabral_spearman <- cor_pmat(corr_Cabral_spearman)

# Specify the desired order of parameters
#param_orderC <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names

# Reorder the correlation matrix
#corr_Cabral2 <- corr_Cabral2[param_orderC, param_orderC]

# Extract the correlations involving "Silicate"
corr_silicate_C <- corr_Cabral_spearman[, "Silicate", drop = FALSE]

# Transpose the matrix to make "Silicate" the row instead of the column
corr_silicate_matrix_C <- t(as.matrix(corr_silicate_C))

# Convert the matrix to a long format
corr_melted_C <- melt(corr_silicate_matrix_C)

# Plot the correlation matrix vertically
Cabral_corr_matrix_spearman <- ggplot(corr_melted_C, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile(color = "white") +  # Adjust width of the tiles to make the plot narrower
  scale_fill_gradient2(low = "steelblue3", high = "tomato3", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Spearman\nCorrelation") +
  geom_text(aes(label = round(value, 2)), color = "black", size = 5) +
  theme_minimal() + 
  theme(axis.text.x = element_text(vjust = 0.5, size = 12, face = "bold", color="black"),
        axis.text.y = element_text(size = 12, face = "bold", color="black")) +  
  coord_flip() +  # Flip coordinates for vertical layout
  labs(x = "", y = "")

Cabral_corr_matrix_spearman

ggsave(here::here("Outputs", "Results", "Final", "Csilicate_corrmatrix.jpg"), 
       width=7, height=12)


########################### 
## patch together 
###########################

#### remove x-axis labels 
p_Ccorr_silicate <- Cabral_corr_matrix_spearman 
p_Vcorr_silicate <- Varari_corr_matrix_spearman 
 
combined_corr_plot_silicate <- p_Ccorr_silicate/p_Vcorr_silicate + plot_layout(guides = 'collect')

combined_corr_plot_silicate + 
  plot_annotation(tag_levels = 'A') &
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"), 
        plot.tag = element_text(size = 18, face = "bold"))

ggsave(here::here("Outputs", "Results", "Final", "CorrMatrices.jpg"), 
       width=10, height=15)

Getting p-values and r-squared coefficients for Spearman correlation plots

#################################
## Varari 
#################################

# Calculate the Spearman correlation and p-values matrix
res_Varari_spearman <- rcorr(as.matrix(corr_matrix_Varari_silicateonly), type = "spearman")

# Extract the correlation matrix
corr_Varari_spearman <- res_Varari_spearman$r

corr_Varari_spearman
##              Salinity          TA         pH        N+N  Phosphate   Silicate
## Salinity   1.00000000  0.06509682 -0.1481712  0.2119188 -0.3863537 -0.3078341
## TA         0.06509682  1.00000000 -0.5024418  0.7087955  0.3498518  0.8119126
## pH        -0.14817116 -0.50244181  1.0000000 -0.5180828 -0.1335543 -0.4445489
## N+N        0.21191882  0.70879554 -0.5180828  1.0000000  0.2393305  0.5343213
## Phosphate -0.38635374  0.34985177 -0.1335543  0.2393305  1.0000000  0.4618148
## Silicate  -0.30783408  0.81191260 -0.4445489  0.5343213  0.4618148  1.0000000
# Extract the p-values matrix
p.mat_Varari_spearman <- res_Varari_spearman$P

p.mat_Varari_spearman
##               Salinity           TA           pH          N+N    Phosphate
## Salinity            NA 5.896440e-01 2.175118e-01 7.603422e-02 8.749647e-04
## TA        0.5896439781           NA 8.015339e-06 4.647394e-12 2.783343e-03
## pH        0.2175117955 8.015339e-06           NA 3.702010e-06 2.668455e-01
## N+N       0.0760342177 4.647394e-12 3.702010e-06           NA 4.441475e-02
## Phosphate 0.0008749647 2.783343e-03 2.668455e-01 4.441475e-02           NA
## Silicate  0.0090129353 0.000000e+00 1.027994e-04 1.592406e-06 5.036852e-05
##               Silicate
## Salinity  9.012935e-03
## TA        0.000000e+00
## pH        1.027994e-04
## N+N       1.592406e-06
## Phosphate 5.036852e-05
## Silicate            NA
# Calculate approximate r-squared by squaring the correlation coefficients
r_squared_Varari_spearman <- corr_Varari_spearman^2

r_squared_Varari_spearman
##              Salinity          TA         pH        N+N  Phosphate   Silicate
## Salinity  1.000000000 0.004237596 0.02195469 0.04490959 0.14926921 0.09476182
## TA        0.004237596 1.000000000 0.25244777 0.50239112 0.12239626 0.65920206
## pH        0.021954692 0.252447770 1.00000000 0.26840979 0.01783675 0.19762370
## N+N       0.044909586 0.502391121 0.26840979 1.00000000 0.05727908 0.28549925
## Phosphate 0.149269214 0.122396261 0.01783675 0.05727908 1.00000000 0.21327287
## Silicate  0.094761820 0.659202062 0.19762370 0.28549925 0.21327287 1.00000000
##########################
## only for silicate

# Extract correlations involving Silicate
corr_silicate_V <- corr_Varari_spearman[, "Silicate", drop = FALSE]

# Extract p-values involving Silicate
p_values_silicate_V <- p.mat_Varari_spearman[, "Silicate", drop = FALSE]

# Extract r-squared values involving Silicate
r_squared_silicate_V <- r_squared_Varari_spearman[, "Silicate", drop = FALSE]

# Print results for Silicate
print(corr_silicate_V)
##             Silicate
## Salinity  -0.3078341
## TA         0.8119126
## pH        -0.4445489
## N+N        0.5343213
## Phosphate  0.4618148
## Silicate   1.0000000
print(p_values_silicate_V)
##               Silicate
## Salinity  9.012935e-03
## TA        0.000000e+00
## pH        1.027994e-04
## N+N       1.592406e-06
## Phosphate 5.036852e-05
## Silicate            NA
print(r_squared_silicate_V)
##             Silicate
## Salinity  0.09476182
## TA        0.65920206
## pH        0.19762370
## N+N       0.28549925
## Phosphate 0.21327287
## Silicate  1.00000000
# Identify significant correlations (p-value < 0.05)
significant_correlations <- p_values_silicate_V < 0.05

# Extract significant correlations
significant_corrs <- corr_silicate_V[significant_correlations]

# Extract the corresponding r-squared values
significant_r_squared <- r_squared_silicate_V[significant_correlations]

# Print significant correlations and their r-squared values
print(significant_corrs)
## [1] -0.3078341  0.8119126 -0.4445489  0.5343213  0.4618148         NA
print(significant_r_squared)
## [1] 0.09476182 0.65920206 0.19762370 0.28549925 0.21327287         NA
#################################
## Cabral 
#################################

# Calculate the Spearman correlation and p-values matrix
res_Cabral_spearman <- rcorr(as.matrix(corr_matrix_Cabral), type = "spearman")

# Extract the correlation matrix
corr_Cabral_spearman <- res_Cabral_spearman$r

corr_Cabral_spearman
##              Salinity          TA          pH         N+N  Phosphate
## Salinity   1.00000000  0.04362344  0.07363389  0.18414095 -0.5645513
## TA         0.04362344  1.00000000  0.02442002 -0.09770993  0.3998131
## pH         0.07363389  0.02442002  1.00000000  0.24030535 -0.1425907
## N+N        0.18414095 -0.09770993  0.24030535  1.00000000 -0.1536658
## Phosphate -0.56455133  0.39981309 -0.14259068 -0.15366585  1.0000000
## Silicate  -0.80452665  0.03785682  0.14257365 -0.18964728  0.5856742
##              Silicate
## Salinity  -0.80452665
## TA         0.03785682
## pH         0.14257365
## N+N       -0.18964728
## Phosphate  0.58567422
## Silicate   1.00000000
# Extract the p-values matrix
p.mat_Cabral_spearman <- res_Cabral_spearman$P

p.mat_Cabral_spearman
##               Salinity          TA         pH        N+N    Phosphate
## Salinity            NA 0.754118246 0.59669822 0.18255726 8.726575e-06
## TA        7.541182e-01          NA 0.86086202 0.48211665 2.741903e-03
## pH        5.966982e-01 0.860862023         NA 0.08006565 3.036804e-01
## N+N       1.825573e-01 0.482116652 0.08006565         NA 2.672614e-01
## Phosphate 8.726575e-06 0.002741903 0.30368044 0.26726141           NA
## Silicate  2.322587e-13 0.785792924 0.30373878 0.16959400 3.283811e-06
##               Silicate
## Salinity  2.322587e-13
## TA        7.857929e-01
## pH        3.037388e-01
## N+N       1.695940e-01
## Phosphate 3.283811e-06
## Silicate            NA
# Calculate approximate r-squared by squaring the correlation coefficients
r_squared_Cabral_spearman <- corr_Cabral_spearman^2

r_squared_Cabral_spearman
##              Salinity           TA           pH        N+N  Phosphate
## Salinity  1.000000000 0.0019030045 0.0054219499 0.03390789 0.31871820
## TA        0.001903004 1.0000000000 0.0005963376 0.00954723 0.15985050
## pH        0.005421950 0.0005963376 1.0000000000 0.05774666 0.02033210
## N+N       0.033907888 0.0095472301 0.0577466635 1.00000000 0.02361319
## Phosphate 0.318718203 0.1598505036 0.0203321023 0.02361319 1.00000000
## Silicate  0.647263135 0.0014331385 0.0203272470 0.03596609 0.34301429
##              Silicate
## Salinity  0.647263135
## TA        0.001433139
## pH        0.020327247
## N+N       0.035966090
## Phosphate 0.343014290
## Silicate  1.000000000
##########################
## only for silicate

# Extract correlations involving Silicate
corr_silicate_C <- corr_Cabral_spearman[, "Silicate", drop = FALSE]

# Extract p-values involving Silicate
p_values_silicate_C <- p.mat_Cabral_spearman[, "Silicate", drop = FALSE]

# Extract r-squared values involving Silicate
r_squared_silicate_C <- r_squared_Cabral_spearman[, "Silicate", drop = FALSE]

# Print results for Silicate
print(corr_silicate_C)
##              Silicate
## Salinity  -0.80452665
## TA         0.03785682
## pH         0.14257365
## N+N       -0.18964728
## Phosphate  0.58567422
## Silicate   1.00000000
print(p_values_silicate_C)
##               Silicate
## Salinity  2.322587e-13
## TA        7.857929e-01
## pH        3.037388e-01
## N+N       1.695940e-01
## Phosphate 3.283811e-06
## Silicate            NA
print(r_squared_silicate_C)
##              Silicate
## Salinity  0.647263135
## TA        0.001433139
## pH        0.020327247
## N+N       0.035966090
## Phosphate 0.343014290
## Silicate  1.000000000
# Identify significant correlations (p-value < 0.05)
significant_correlations <- p.mat_Cabral_spearman < 0.05

# Extract significant correlations
significant_corrs <- corr_Cabral_spearman[significant_correlations]

# Extract the corresponding r-squared values
significant_r_squared <- r_squared_Cabral_spearman[significant_correlations]

# Print significant correlations and their r-squared values
print(significant_corrs)
##  [1]         NA -0.5645513 -0.8045267         NA  0.3998131         NA
##  [7]         NA -0.5645513  0.3998131         NA  0.5856742 -0.8045267
## [13]  0.5856742         NA
print(significant_r_squared)
##  [1]        NA 0.3187182 0.6472631        NA 0.1598505        NA        NA
##  [8] 0.3187182 0.1598505        NA 0.3430143 0.6472631 0.3430143        NA

Figure 2a (part2): Barplots for GW endmember data

  • making bar graphs to illustrate differences between sites and from sites to offshore ambient SW
###########################################################
### making bar graphs to illustrate differences ### 
##########################################################
 
#### SEE BACKGROUND SITE DATA SCRIPT FOR HOW OFFSHORE VALUES WERE CALCULATED ####

############################## 
### plots for CSUNposium 
###############################
sitecomparison_maxdata <- tribble(~Site, ~Parameter, ~Value, 
        "Varari", "pH", "7.73", 
        "Cabral", "pH", "6.89", 
        "Reef Crest", "pH", "7.9", 
        "Varari", "Nitrates + Nitrites (umol/L)", "42.49", 
        "Cabral", "Nitrates + Nitrites (umol/L)", "8.22", 
        "Reef Crest", "Nitrates + Nitrites (umol/L)", "0.58", 
        "Varari", "TA (umol/kg-1)", "5393.3", 
        "Cabral", "TA (umol/kg-1)", "1755.18", 
        "Reef Crest", "TA (umol/kg-1)", "2386.84", 
        "Varari", "Phosphate (umol/L)", "2.5", 
        "Cabral", "Phosphate (umol/L)", "3.61", 
        "Reef Crest", "Phosphate (umol/L)", "0.20", 
        "Varari", "Silicate (umol/L)", "725.61", 
        "Cabral", "Silicate (umol/L)", "855.01", 
        "Reef Crest", "Silicate (umol/L)", "3.68", 
        "Varari", "Salinity (psu)", "6.32", 
        "Cabral", "Salinity (psu)", "5.79", 
        "Reef Crest", "Salinity (psu)", "36.7") %>% 
  mutate(Site= as.factor(Site))


### for pH 
SiteComparisons_pH <- sitecomparison_maxdata %>% 
  filter(Parameter=="pH") %>% 
  ggplot(aes(x=Site,
             y=Value, 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  labs(x="Location", 
       y = expression(pH[T])) +
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_pH

ggsave(here::here("Outputs", "Results", "SitevOffshore", "pH.jpg"))

## for TA 
SiteComparisons_TA <- sitecomparison_maxdata %>% 
  filter(Parameter=="TA (umol/kg-1)")  %>% 
  ggplot(aes(x=Site,
             y=Value, 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  labs(x="Location",
       y = expression(TA~(mu*mol~kg^-1))) +
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_TA

ggsave(here::here("Outputs", "Results", "SitevOffshore", "TA.jpg"))

## Salinity 
SiteComparisons_Salinity <- sitecomparison_maxdata %>% 
  filter(Parameter=="Salinity (psu)") %>% 
  ggplot(aes(x=Site,
             y=Value, 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_y_discrete(limits = c("5.79", "6.32", "36.7")) +
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  labs(x="Location", 
       y = "Salinity (psu)") +
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_Salinity

ggsave(here::here("Outputs", "Results", "SitevOffshore", "Salinity.jpg"))

## for Phosphates 
SiteComparisons_Phosphate <- sitecomparison_maxdata %>% 
  filter(Parameter=="Phosphate (umol/L)")  %>% 
  ggplot(aes(x=Site,
             y=Value, 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  labs(x="Location", 
       y = expression(Phosphate~(mu*mol~L^-1))) +
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_Phosphate

ggsave(here::here("Outputs", "Results", "SitevOffshore", "Phosphate.jpg"))

## for Silicate 
SiteComparisons_Silicate <- sitecomparison_maxdata %>% 
  filter(Parameter=="Silicate (umol/L)")  %>% 
  ggplot(aes(x=Site,
             y=Value, 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  labs(x="Location", 
       y = expression(Silicate~(mu*mol~L^-1))) + 
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_Silicate

ggsave(here::here("Outputs", "Results", "SitevOffshore", "Silicate.jpg"))

  
## nitrates 
SiteComparisons_Nitrates <- sitecomparison_maxdata %>% 
  filter(Parameter=="Nitrates + Nitrites (umol/L)") %>% 
  ggplot(aes(x=Site,
             y=as.numeric(Value), 
             fill=Site)) + 
  geom_col()  + 
  scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) + 
  scale_y_continuous(trans='log10') +
  #scale_y_continuous(trans = "log10", # Apply logarithmic transformation
           #         breaks = c(0.58, 8.22, 42.49),  # Set specific breaks
              #     labels = c("0.58", "8.22", "42.49")) +  # Custom labels for clarity
  scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"), 
                    limits=c("Reef Crest", "Cabral", "Varari"), 
                    guide="none") +
  theme_classic() + 
  #scale_y_continuous(trans='log10') +
  labs(y = expression(Nitrate+Nitrite~(mu*mol~L^-1))) +
  theme(axis.text.x=element_text(size=25), 
        axis.text.y=element_text(size=25), 
        axis.title.x=element_text(size=25), 
        axis.title.y=element_text(size=25))

SiteComparisons_Nitrates

ggsave(here::here("Outputs", "Results", "SitevOffshore", "NN.jpg"))


### patch them all together 
p1.1 <- SiteComparisons_Nitrates + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.2 <- SiteComparisons_Phosphate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.3 <- SiteComparisons_Silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.4 <- SiteComparisons_Salinity + theme(axis.title.x = element_blank()) 
p1.5<- SiteComparisons_TA + theme(axis.title.x = element_blank()) 
p1.6 <- SiteComparisons_pH + theme(axis.title.x = element_blank()) 

combined_endmemberSW_plot <- (p1.1 + p1.2 + p1.3) / (p1.4 + p1.5 + p1.6) +
  plot_layout(guides = 'collect') + 
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 24, face = "bold"))

combined_endmemberSW_plot

ggsave(here::here("Outputs", "Results", "Final", "rawGW_barplots.jpg"), 
       width=20, height = 12)



#plot.margin = unit(c(1, 1, 1, 1), "cm"), 

Figure 2b: Background information (plots with silicate vs parameters)

  • for personal understanding - not to include in paper
##########################
## site data  
##########################
background_site_plots <- RespoR_Normalized_Full_withNEC_nutrients %>% 
 # mutate(N_P = NN_umolL/phosphate_umolL) %>% 
  filter(P_R=="C") %>% 
  mutate("Salinity"=salinity, 
         "TA" = TA_initial, 
         "pH" = new_pH, 
         "NN" = NN_umolL, 
         "Phosphate" = phosphate_umolL,
         "Silicate" = silicate_umolL) %>% 
  select(!c("salinity", "log_silicate", "log_silicate2", "silicate_umolL", "TA_initial", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "date", "P_R", "ammonia_umolL")) %>% 
  filter(Phosphate<0.9 & TA < 2550)

##########################
## ggplots - silciate vs everything else 
##########################

### Phos 
silicate_vs_phos <- background_site_plots %>% 
  ggplot(aes(x=Silicate, 
             y=Phosphate, 
             color=site)) + 
  geom_point(size=3) + 
  facet_wrap(~site, scales = "free_y") ## OUTLIER AT 0.9 for PHOS (removed)

silicate_vs_phos

### NN
silicate_vs_NN <- background_site_plots %>% 
  ggplot(aes(x=Silicate, 
             y= NN,
             color=site)) + 
  geom_point(size=3) + 
  facet_wrap(~site, scales = "free_y") 

silicate_vs_NN

### TA
silicate_vs_TA <- background_site_plots %>% 
  ggplot(aes(x= Silicate, 
             y= TA,
             color=site)) + 
  geom_point(size=3) + 
  facet_wrap(~site, scales = "free_y") ## OUTLIER AT VARARI ABOVE 2550 

silicate_vs_TA

### Salinity
silicate_vs_salinity <- background_site_plots %>% 
  ggplot(aes(x=Silicate, 
             y= Salinity,
             color=site)) + 
  geom_point(size=3) + 
  facet_wrap(~site, scales = "free_y") 

silicate_vs_salinity

### pH
silicate_vs_pH <- background_site_plots %>% 
  ggplot(aes(x=Silicate, 
             y= pH,
             color=site)) + 
  geom_point(size=3) + 
  facet_wrap(~site, scales = "free_y") 

silicate_vs_pH

############ patch together 
silicate_by_otherparams <- silicate_vs_pH + silicate_vs_salinity + silicate_vs_TA + silicate_vs_NN + silicate_vs_phos

silicate_by_otherparams

ggsave(here::here("Outputs", "Results", "Final", "silicate_vs_otherparams.jpg"), 
       width=24, height = 12)

Figure 3: Physio rate plots

### Respiration 
R_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=site)) +
  scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "R", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y~poly(log(x),2), color = "black") +
  labs(y = expression(Respiration~Rate~(mu*mol~O[2]~cm^-2~hr^-1)),
       color = "Site") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=15, face="bold"),
        axis.title.y=element_text(size=15, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
R_silicate

ggsave(here::here("Outputs", "Results", "Final", "R_silicate.jpg"), 
       width=10, height=7)

RespoR_Normalized_Full_withNEC_nutrientsR <- RespoR_Normalized_Full_withNEC_nutrients %>% 
   filter(P_R=="C") %>% 
 filter(site=="Cabral") %>% 
  filter(!is.infinite(umol.cm2.hr))
  
max(RespoR_Normalized_Full_withNEC_nutrientsR$umol.cm2.hr)
## [1] 0.4004983
min(RespoR_Normalized_Full_withNEC_nutrientsR$umol.cm2.hr)
## [1] -0.06829634
###############################
## Gross Photosynthesis 
###############################
GP_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=site)) +
 scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "GP", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y~poly(log(x),2), color = "black") +
  labs(x = expression(Silicate~(mu*mol~L^-1)),
      y = expression(Gross~Photosynthesis~Rate~(mu*mol~O[2]~cm^-2~hr^-1)),
       color = "Site") +
    theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20, face="bold"),
        axis.title.y=element_text(size=20, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
GP_silicate

ggsave(here::here("Outputs", "Results", "GP_silicate.jpg"), 
       width=10, height=7)

###############################
## Calcification 
###############################
C_silicate_ECRS_firebrick4 <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=5, aes(color=site)) +
  scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  #scale_x_discrete(breaks = c(1, 10, 20)) +
  geom_hline(yintercept = 0) +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "C", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y~log(x), color = "black") +
 labs(x = expression(Silicate~(mu*mol~L^-1)),
       y = expression(Calcification~Rate~(mu*mol~CaCO[3]~cm^-2~hr^-1)),
       color = "Site") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=15, face="bold"),
        axis.title.y=element_text(size=15, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
C_silicate_ECRS_firebrick4

ggsave(here::here("Outputs", "Results", "C_silicate.jpg"), 
       width=10, height=7)

#############################################
### make these plots figure worthy using patchwork 
#############################################

#### remove x-axis labels 
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),  
                         legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
                         axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p2 <- GP_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), 
                          axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p3 <- C_silicate_ECRS_firebrick4 + 
  theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
        axis.title.x=element_text(size=15, face="bold"), 
        axis.title.y=element_text(size=15, face="bold")) +
  scale_x_continuous(breaks = c(2, 10, 20)) 

combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')

combined_column + 
  plot_annotation(tag_levels = 'A') &
 theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        plot.tag = element_text(size = 24, face = "bold"))

ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"), 
       width=10, height=16)

Figure 3: Re’Vampd with GP and R nonlinear

###################
### GP 2.0
###################

#### confirm models are significant for both sites 

### confirm nonlinear model 
GP_Varari_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  #filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                         Sum Sq  Mean Sq NumDF   DenDF F value
## poly(log(silicate_umolL + 1), 2)      0.210557 0.105279     2 106.556  6.2548
## site                                  0.009422 0.009422     1  11.731  0.5598
## poly(log(silicate_umolL + 1), 2):site 0.045241 0.022621     2 106.556  1.3439
##                                         Pr(>F)   
## poly(log(silicate_umolL + 1), 2)      0.002701 **
## site                                  0.469077   
## poly(log(silicate_umolL + 1), 2):site 0.265202   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(GP_Varari_silicate_nonlinearmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ poly(log(silicate_umolL + 1), 2) * site + (1 |  
##     new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -120.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9126 -0.6755  0.1116  0.6361  3.4111 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  new_colonynumber (Intercept) 0.00760  0.08718 
##  Residual                     0.01683  0.12974 
## Number of obs: 122, groups:  new_colonynumber, 14
## 
## Fixed effects:
##                                               Estimate Std. Error        df
## (Intercept)                                    0.76218    0.03994  11.72728
## poly(log(silicate_umolL + 1), 2)1             -0.53415    0.22012 104.47568
## poly(log(silicate_umolL + 1), 2)2              0.43175    0.23982 105.71662
## siteVarari                                    -0.03954    0.05285  11.73145
## poly(log(silicate_umolL + 1), 2)1:siteVarari   0.43517    0.28499 105.40538
## poly(log(silicate_umolL + 1), 2)2:siteVarari   0.08519    0.30228 108.77719
##                                              t value Pr(>|t|)    
## (Intercept)                                   19.081 3.42e-10 ***
## poly(log(silicate_umolL + 1), 2)1             -2.427   0.0169 *  
## poly(log(silicate_umolL + 1), 2)2              1.800   0.0747 .  
## siteVarari                                    -0.748   0.4691    
## poly(log(silicate_umolL + 1), 2)1:siteVarari   1.527   0.1298    
## poly(log(silicate_umolL + 1), 2)2:siteVarari   0.282   0.7786    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##               (Intr) pl((_L+1),2)1 pl((_L+1),2)2 sitVrr p((_L+1),2)1:
## pl((_L+1),2)1 -0.068                                                 
## pl((_L+1),2)2  0.070 -0.413                                          
## siteVarari    -0.756  0.052        -0.053                            
## p((_L+1),2)1:  0.053 -0.772         0.319        -0.026              
## p((_L+1),2)2: -0.056  0.328        -0.793         0.029 -0.199
r.squaredGLMM(GP_Varari_silicate_nonlinearmodel)
##             R2m       R2c
## [1,] 0.09963739 0.3797085
#### now plot 
GP_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=site)) +
 scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "GP", site == "Varari", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y~poly(log(x),2), color = "firebrick4") +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "GP", site == "Cabral", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y~poly(log(x),2), color = "darkgoldenrod2") +
  labs(x = expression(Silicate~(mu*mol~L^-1)),
      y = expression(Gross~Photosynthesis~Rate~(mu*mol~O[2]~cm^2~hr^-1)),
       color = "Site") +
    theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20, face="bold"),
        axis.title.y=element_text(size=20, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
GP_silicate2 

########################### 
## patch together 
###########################

#### remove x-axis labels 
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),  
                         legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
                         axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p2 <- GP_silicate2 + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), 
                          axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p3 <- C_silicate_ECRS_firebrick4 + 
  theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
        axis.title.x=element_text(size=15, face="bold"), 
        axis.title.y=element_text(size=15, face="bold")) +
  scale_x_continuous(breaks = c(2, 10, 20)) 

combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')

combined_column + 
  plot_annotation(tag_levels = 'A') &
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        plot.tag = element_text(size = 18, face = "bold"))

ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"), 
       width=10, height=15)

FINAL MODELS USED IN RESULTS:

RESULTS SHOW THAT: - linear model is better fit for R, GP, and C - decided on ANCOVA design rather than analyzing sites separately – acknowledging that they are different BUT want to test for the relationship between the two to see HOW they are different and how the GP, R, and C rates are different in the fragments exposed to gw from both sites - model set-up: “model <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=data” – and then filter for specifically R, GP, or C

  • Respiration: log_silicate, 2 (poly) = 0.0065 and log_silicate is 0.0045, so both significant. The interaction between silicate and site = 0.95 and just site (0.2) are both insignificant
  • Photosynthesis: Poly (0.001) and linear silicate (0.0005) both significant. Interaction and just site are both insignificant
  • Calcification: ONLY testing linear. Logged silicate = 0.01 so significant. Interaction and just site are both insignificant
RespoR_Normalized_Full_withNEC_nutrients <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  mutate(log_silicate = log(silicate_umolL+1), 
         log_silicate2 = (log_silicate^2), 
         date=factor(date))
########################
## R
########################
### nonlinear model 
#R_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
 # filter(P_R=="R") %>% 
  #filter(!is.infinite(umol.cm2.hr)))

R_silicate_nonlinearmodel2 <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_silicate_nonlinearmodel2)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)   
## log_silicate      0.046571 0.046571     1 107.262  6.9107 0.009826 **
## site              0.000212 0.000212     1  47.521  0.0315 0.859918   
## log_silicate:site 0.002363 0.002363     1 107.262  0.3507 0.554987   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(R_silicate_nonlinearmodel2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + (1 | new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "C") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -214.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6824 -0.5040 -0.0080  0.5366  2.8061 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  new_colonynumber (Intercept) 0.006216 0.07884 
##  Residual                     0.006739 0.08209 
## Number of obs: 122, groups:  new_colonynumber, 14
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)               0.16360    0.05017  48.53078   3.261  0.00204 **
## log_silicate             -0.02525    0.01844 106.08033  -1.369  0.17385   
## siteVarari               -0.01173    0.06611  47.52070  -0.177  0.85992   
## log_silicate:siteVarari  -0.01468    0.02480 107.26152  -0.592  0.55499   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr
## log_silicat -0.733              
## siteVarari  -0.759  0.557       
## lg_slct:stV  0.545 -0.744 -0.730
check_model(R_silicate_nonlinearmodel2)

AIC(R_silicate_nonlinearmodel2)
## [1] -202.3236
r.squaredGLMM(R_silicate_nonlinearmodel2)
##             R2m       R2c
## [1,] 0.05804063 0.5099971
### linear model 
#R_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  #filter(site=="Varari") %>% 
 # filter(P_R=="R") %>% 
 # filter(!is.infinite(umol.cm2.hr)))

R_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
                                            filter(P_R=="R") %>% 
                                            filter(!is.infinite(umol.cm2.hr)))
anova(R_bothsites_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq   Mean Sq NumDF   DenDF F value Pr(>F)
## log_silicate      0.0042170 0.0042170     1 107.372  1.2479 0.2664
## site              0.0018190 0.0018190     1  53.797  0.5383 0.4663
## log_silicate:site 0.0014839 0.0014839     1 107.372  0.4391 0.5090
summary(R_bothsites_silicate_linearmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + (1 | new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -297.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.77505 -0.58754  0.08962  0.56703  2.36504 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  new_colonynumber (Intercept) 0.002584 0.05083 
##  Residual                     0.003379 0.05813 
## Number of obs: 122, groups:  new_colonynumber, 14
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)               0.361900   0.034255  55.124914  10.565  7.3e-15 ***
## log_silicate             -0.003988   0.013059 105.976258  -0.305    0.761    
## siteVarari               -0.033095   0.045107  53.796820  -0.734    0.466    
## log_silicate:siteVarari  -0.011629   0.017549 107.372117  -0.663    0.509    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr
## log_silicat -0.761              
## siteVarari  -0.759  0.578       
## lg_slct:stV  0.566 -0.744 -0.757
AIC(R_bothsites_silicate_linearmodel)
## [1] -285.7517
########################
## GP
########################
### nonlinear model 
GP_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)    
## log_silicate      0.212785 0.212785     1 112.396 12.7387 0.000528 ***
## site              0.051338 0.051338     1  80.138  3.0734 0.083404 .  
## log_silicate2     0.185984 0.185984     1 111.720 11.1343 0.001151 ** 
## log_silicate:site 0.043913 0.043913     1 105.552  2.6289 0.107916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(GP_silicate_nonlinearmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## umol.cm2.hr ~ log_silicate * site + log_silicate2 + (1 | new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -110.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9533 -0.6610  0.1187  0.6322  3.4536 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  new_colonynumber (Intercept) 0.007516 0.0867  
##  Residual                     0.016704 0.1292  
## Number of obs: 122, groups:  new_colonynumber, 14
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               1.39873    0.17566 116.51375   7.963 1.24e-12 ***
## log_silicate             -0.57275    0.15819 111.57997  -3.621 0.000444 ***
## siteVarari               -0.16698    0.09525  80.13835  -1.753 0.083404 .  
## log_silicate2             0.11436    0.03427 111.71969   3.337 0.001151 ** 
## log_silicate:siteVarari   0.06559    0.04045 105.55218   1.621 0.107916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr lg_sl2
## log_silicat -0.962                     
## siteVarari  -0.524  0.358              
## log_silict2  0.917 -0.983 -0.250       
## lg_slct:stV  0.485 -0.398 -0.834  0.271
AIC(GP_silicate_nonlinearmodel)
## [1] -96.5832
r.squaredGLMM(GP_silicate_nonlinearmodel)
##             R2m       R2c
## [1,] 0.09906217 0.3786475
### linear model 
#GP_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  #filter(site=="Varari") %>% 
 # filter(P_R=="GP") %>% 
  #filter(!is.infinite(umol.cm2.hr)))

#anova(GP_silicate_linearmodel)
#summary(GP_silicate_linearmodel)
#AIC(GP_silicate_linearmodel)

########################
## Calcification 
########################
### nonlinear model 
#C_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  #filter(site=="Varari") %>% 
#  filter(P_R=="C") %>% 
#  filter(!is.infinite(umol.cm2.hr)))

#anova(C_Varari_silicate_nonlinearmodel)
#summary(C_Varari_silicate_nonlinearmodel)
#AIC(C_Varari_silicate_nonlinearmodel)

### linear model 
C_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)   
## log_silicate      0.046571 0.046571     1 107.262  6.9107 0.009826 **
## site              0.000212 0.000212     1  47.521  0.0315 0.859918   
## log_silicate:site 0.002363 0.002363     1 107.262  0.3507 0.554987   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(C_Varari_silicate_linearmodel)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log(silicate_umolL + 1) * site + (1 | new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "C") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -214.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6824 -0.5040 -0.0080  0.5366  2.8061 
## 
## Random effects:
##  Groups           Name        Variance Std.Dev.
##  new_colonynumber (Intercept) 0.006216 0.07884 
##  Residual                     0.006739 0.08209 
## Number of obs: 122, groups:  new_colonynumber, 14
## 
## Fixed effects:
##                                     Estimate Std. Error        df t value
## (Intercept)                          0.16360    0.05017  48.53078   3.261
## log(silicate_umolL + 1)             -0.02525    0.01844 106.08033  -1.369
## siteVarari                          -0.01173    0.06611  47.52070  -0.177
## log(silicate_umolL + 1):siteVarari  -0.01468    0.02480 107.26152  -0.592
##                                    Pr(>|t|)   
## (Intercept)                         0.00204 **
## log(silicate_umolL + 1)             0.17385   
## siteVarari                          0.85992   
## log(silicate_umolL + 1):siteVarari  0.55499   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg(_L+1) sitVrr
## lg(slc_L+1) -0.733                
## siteVarari  -0.759  0.557         
## lg(s_L+1):V  0.545 -0.744   -0.730
AIC(C_Varari_silicate_linearmodel)
## [1] -202.3236
r.squaredGLMM(C_silicate_linearmodel)
##             R2m       R2c
## [1,] 0.05804063 0.5099971
########################################################
########################################################

Models worked on August 29th, 2024

#####################
## Respiration 
#####################
modeltest_R <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="R") %>%
  filter(!is.infinite(umol.cm2.hr)))

anova(modeltest_R)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq   Mean Sq NumDF   DenDF F value   Pr(>F)   
## log_silicate      0.0290541 0.0290541     1 106.795  9.2865 0.002909 **
## site              0.0020049 0.0020049     1   9.596  0.6408 0.442786   
## I(log_silicate^2) 0.0265717 0.0265717     1 106.625  8.4931 0.004345 **
## log_silicate:site 0.0000016 0.0000016     1 105.325  0.0005 0.981786   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modeltest_R)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |  
##     date/new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -307.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.49031 -0.60178  0.00364  0.73007  2.41943 
## 
## Random effects:
##  Groups                Name        Variance  Std.Dev.
##  new_colonynumber:date (Intercept) 0.0001486 0.01219 
##  date                  (Intercept) 0.0044872 0.06699 
##  Residual                          0.0031286 0.05593 
## Number of obs: 122, groups:  new_colonynumber:date, 14; date, 7
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)              5.683e-01  8.496e-02  6.160e+01   6.690 7.59e-09 ***
## log_silicate            -2.032e-01  6.949e-02  1.066e+02  -2.924  0.00422 ** 
## siteVarari              -5.035e-02  6.290e-02  9.596e+00  -0.801  0.44279    
## I(log_silicate^2)        4.389e-02  1.506e-02  1.066e+02   2.914  0.00434 ** 
## log_silicate:siteVarari  4.011e-04  1.753e-02  1.053e+02   0.023  0.98179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.873                     
## siteVarari  -0.548  0.231              
## I(lg_slc^2)  0.833 -0.984 -0.162       
## lg_slct:stV  0.428 -0.385 -0.545  0.260
#####################
## Photosynthesis  
#####################
modeltest_GP <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="GP") %>%
  filter(!is.infinite(umol.cm2.hr)))

anova(modeltest_GP)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq  Mean Sq NumDF   DenDF F value    Pr(>F)    
## log_silicate      0.220744 0.220744     1 109.495 13.2207 0.0004236 ***
## site              0.050153 0.050153     1  36.857  3.0037 0.0914350 .  
## I(log_silicate^2) 0.193149 0.193149     1 109.259 11.5680 0.0009381 ***
## log_silicate:site 0.042678 0.042678     1 105.478  2.5560 0.1128650    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modeltest_GP)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |  
##     date/new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -110.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9118 -0.6580  0.1101  0.6406  3.4806 
## 
## Random effects:
##  Groups                Name        Variance Std.Dev.
##  new_colonynumber:date (Intercept) 0.005917 0.07692 
##  date                  (Intercept) 0.001941 0.04406 
##  Residual                          0.016697 0.12922 
## Number of obs: 122, groups:  new_colonynumber:date, 14; date, 7
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               1.41148    0.17736 107.13408   7.958 1.94e-12 ***
## log_silicate             -0.58502    0.15884 109.21295  -3.683 0.000360 ***
## siteVarari               -0.17129    0.09883  36.85738  -1.733 0.091435 .  
## I(log_silicate^2)         0.11705    0.03442 109.25888   3.401 0.000938 ***
## log_silicate:siteVarari   0.06469    0.04046 105.47770   1.599 0.112865    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.956                     
## siteVarari  -0.523  0.345              
## I(lg_slc^2)  0.912 -0.983 -0.242       
## lg_slct:stV  0.478 -0.394 -0.803  0.268

PLOTS BASED ON MODELS: RESPIRATION

###################
### R 
###################

modeltest_R <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="R") %>%
  filter(!is.infinite(umol.cm2.hr)))

anova(modeltest_R)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq   Mean Sq NumDF   DenDF F value   Pr(>F)   
## log_silicate      0.0290541 0.0290541     1 106.795  9.2865 0.002909 **
## site              0.0020049 0.0020049     1   9.596  0.6408 0.442786   
## I(log_silicate^2) 0.0265717 0.0265717     1 106.625  8.4931 0.004345 **
## log_silicate:site 0.0000016 0.0000016     1 105.325  0.0005 0.981786   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modeltest_R)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |  
##     date/new_colonynumber)
##    Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%  
##     filter(!is.infinite(umol.cm2.hr))
## 
## REML criterion at convergence: -307.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.49031 -0.60178  0.00364  0.73007  2.41943 
## 
## Random effects:
##  Groups                Name        Variance  Std.Dev.
##  new_colonynumber:date (Intercept) 0.0001486 0.01219 
##  date                  (Intercept) 0.0044872 0.06699 
##  Residual                          0.0031286 0.05593 
## Number of obs: 122, groups:  new_colonynumber:date, 14; date, 7
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)              5.683e-01  8.496e-02  6.160e+01   6.690 7.59e-09 ***
## log_silicate            -2.032e-01  6.949e-02  1.066e+02  -2.924  0.00422 ** 
## siteVarari              -5.035e-02  6.290e-02  9.596e+00  -0.801  0.44279    
## I(log_silicate^2)        4.389e-02  1.506e-02  1.066e+02   2.914  0.00434 ** 
## log_silicate:siteVarari  4.011e-04  1.753e-02  1.053e+02   0.023  0.98179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.873                     
## siteVarari  -0.548  0.231              
## I(lg_slc^2)  0.833 -0.984 -0.162       
## lg_slct:stV  0.428 -0.385 -0.545  0.260
### based on reference from stack overflow on how to visualize a mixed model with random effects 
### create a new dataframe that has x var (silicate), z var (color by colony number), and then create y values using predict function
### logged or not logged silicate? 

### first isolate just R values 
R_only <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  filter(P_R=="R") 

write_csv(R_only, here::here("Data", "R_only.csv"))
  

### now make new df 
plot_df_R <- data.frame(silicate = R_only$silicate_umolL,
                        pred_umol = predict(modeltest_R, newdata = R_only), 
                        new_colonynumber = as.factor(R_only$new_colonynumber))


### define the desired color groups based on colony number (z variable)
plot_df_R$color_group <- ifelse(plot_df_R$new_colonynumber %in% as.character(1:8), "red2", 
                                ifelse(plot_df_R$new_colonynumber %in% as.character(13:18), "gold2", NA))

### define the gradients (Varari in red, Cabral in blue)
Varari_gradient <- scales::gradient_n_pal(c("red2", "coral4"))(seq(0, 1, length.out = length(unique(plot_df_R$new_colonynumber[plot_df_R$color_group == "red2"]))))

Cabral_gradient <- scales::gradient_n_pal(c("gold2", "goldenrod"))(seq(0, 1, length.out = length(unique(plot_df_R$new_colonynumber[plot_df_R$color_group == "gold2"]))))

### Create the ggplot with conditional gradients
Rplot <- plot_df_R %>% 
  ggplot(aes(x = silicate, 
             y = pred_umol,
             color = new_colonynumber)) +
  geom_point(size = 3, alpha = 0.5, 
             data = R_only, 
             aes(x= silicate_umolL,
                 y= umol.cm2.hr, 
                 color = as.factor(new_colonynumber))) +
  geom_line(data = plot_df_R, linetype = 1) +
  coord_trans(x ="log") +
 # geom_abline(aes(intercept = fixef(modeltest_R)[1], slope = fixef(modeltest_R)[2],
           #       linetype = "Fixed effect"), key_glyph = draw_key_path,
           #   linewidth = 1) +
  scale_color_manual(values = c(Varari_gradient, Cabral_gradient)) +
  labs(linetype = NULL, color = "Colony Number", color_group = "Site") +
 # annotation_custom(grid::textGrob(
    #bquote(bold(~y == 
          #   .(round(fixef(modeltest_R)[1], 3)) + .(round(fixef(modeltest_R)[2], 3)) * x)),
   # x = 0.05, y = 0.9, hjust = 0, gp = grid::gpar(cex = 1.2))) +
  ggtitle("Illustration of mixed effects model with random intercept for RESPIRATION") +
  theme(legend.position = "top",
        plot.title = element_text(face = 2, hjust = 0.5),
        panel.border = element_rect(fill = NA, linewidth = 0.2))

Rplot

PLOTS BASED ON MODELS: PHOTOSYNTHESIS

###################
### GP
###################

### based on reference from stack overflow on how to visualize a mixed model with random effects 
### create a new dataframe that has x var (silicate), z var (color by colony number), and then create y values using predict function
### logged or not logged silicate? 

### first isolate just R values 
GP_only <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  filter(P_R=="GP")  

### now make new df 
plot_df_GP <- data.frame(silicate = GP_only$log_silicate,
                        pred_umol = predict(modeltest_GP, newdata = GP_only), 
                        new_colonynumber = as.factor(GP_only$new_colonynumber))


### define the desired color groups based on colony number (z variable)
plot_df_GP$color_group <- ifelse(plot_df_GP$new_colonynumber %in% as.character(1:8), "firebrick4", 
                              ifelse(plot_df_GP$new_colonynumber %in% as.character(13:18), "goldenrod", NA))

### define the gradients (Varari in red, Cabral in blue)
red_gradient <- scales::gradient_n_pal(c("firebrick4"))(seq(0, 1, length.out = length(unique(plot_df_GP$new_colonynumber[plot_df_GP$color_group == "firebrick4"]))))

blue_gradient <- scales::gradient_n_pal(c("goldenrod"))(seq(0, 1, length.out = length(unique(plot_df_GP$new_colonynumber[plot_df_GP$color_group == "goldenrod"]))))

### Create the ggplot with conditional gradients
GPplot <- plot_df_GP %>% 
  ggplot(aes(x = silicate, 
             y = pred_umol,
             color = new_colonynumber)) +
  geom_point(size = 3, alpha = 0.5, 
             data = GP_only, 
             aes(x= log_silicate,
                 y= umol.cm2.hr, 
                 color = as.factor(new_colonynumber))) +
  geom_line(data = plot_df_GP, linetype = 2) +
 # geom_abline(aes(intercept = fixef(modeltest_GP)[1], slope = fixef(modeltest_GP)[2],
              #    linetype = "Fixed effect"), key_glyph = draw_key_path,
             # linewidth = 1) +
  scale_color_manual(values = c(red_gradient, blue_gradient)) +
  labs(linetype = NULL, color = "Colony Number") +
  annotation_custom(grid::textGrob(
    bquote(bold(~y == 
             .(round(fixef(modeltest_GP)[1], 3)) + .(round(fixef(modeltest_GP)[2], 3)) * x)),
    x = 0.05, y = 0.9, hjust = 0, gp = grid::gpar(cex = 1.2))) +
  ggtitle("Illustration of mixed effects model with random intercept for PHOTOSYNTHESIS") +
  theme(legend.position = "top",
        plot.title = element_text(face = 2, hjust = 0.5),
        panel.border = element_rect(fill = NA, linewidth = 0.2))

GPplot

PLOTS BASED ON MODELS: CALCIFICATION

C_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  filter(SGD_number!=0) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr, 
             color=as.factor(new_colonynumber))) + 
  geom_point(size=4) +
  facet_wrap(~new_colonynumber) +
 # scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  geom_smooth(method="lm", formula="y~x + I(x^2)") + 
  labs(x = expression(Silicate~(mu*mol~L^-1)),
      y = expression(Calcification~Rate~(mu*mol~O[2]~cm^2~hr^-1))) +
    theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20, face="bold"),
        axis.title.y=element_text(size=20, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
C_silicate2 

Patch plots together

########################### 
## patch together 
###########################

#### remove x-axis labels 
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),  
                         legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
                         axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p2 <- GP_silicate2 + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), 
                          axis.title.y=element_text(size=15, face="bold")) + 
  scale_x_continuous(breaks = c(2, 10, 20)) 
p3 <- C_silicate_ECRS_firebrick4 + 
  theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none", 
        axis.title.x=element_text(size=15, face="bold"), 
        axis.title.y=element_text(size=15, face="bold")) +
  scale_x_continuous(breaks = c(2, 10, 20)) 

combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')

combined_column + 
  plot_annotation(tag_levels = 'A') &
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        plot.tag = element_text(size = 18, face = "bold"))

ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"), 
       width=10, height=15)

Plots using pred() function

  • as of August 2024, still need to figure out how to do properly
  • as of May 15th, 2025, ignore because did not use in paper
#####################################
## New attempt: July 29th, 2024
#####################################

### new dfs just for calcification 
RespoR_Normalized_Full_withNEC_nutrients_Conly <- RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(P_R=="C") %>%
  filter(!is.infinite(umol.cm2.hr))

### and sort in the right order from lowest silicate to highest 
RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted <- RespoR_Normalized_Full_withNEC_nutrients_Conly %>% 
  arrange(silicate_umolL)

### this is the MIXED MODEL for C
C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), 
                                         data=RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted)

###############################################
## the problem is with the predict function 
###############################################

### try with Riley's code 

new_data <- expand.grid(
  silicate_umolL = seq(min(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL),
                       max(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL), length.out = 100),
  site = levels(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$site),
  new_colonynumber = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$new_colonynumber))

# Predict values using the model
#new_data$predicted <- predict(C_bothsites_silicate_linearmodel, newdata = new_data, re.form = NA)

riley_pred_plot <- new_data %>% 
  ggplot(aes(x = silicate_umolL, 
             y = predicted, 
             color = site)) + 
  geom_line(size = 1, linetype = "dashed") +  # Add dashed trend lines
  geom_point(data = RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted, 
             aes(x = silicate_umolL, 
                 y = umol.cm2.hr, 
                 color = treatment), 
             size = 3) +  # Add original data points
  labs(x = "silicate",
       y = "rate",
       color = "site") +
  theme_minimal() +
  theme(legend.title = element_text(face = "bold"))

#riley_pred_plot

###############################
## tidyverse code
###############################

# Generate a sequence of silicate_umolL values for prediction
new_data_C_pred <- tibble(
    silicate_umolL = seq(min(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL), 
                         max(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL), length.out = 100),
    site = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$site)[2],  # Use the first site as a placeholder
    new_colonynumber = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$new_colonynumber)[5]  # Use the first colony as a placeholder
)

# Add predictions to the new data frame
new_data_C_pred <- new_data_C_pred %>%
    mutate(predicted = predict(C_bothsites_silicate_linearmodel, newdata = new_data_C_pred, re.form = NA))  # predict without random effects


Pred_plot_C_10 <- RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted %>% 
  ggplot(aes(x = log(silicate_umolL), 
             y = umol.cm2.hr, 
             color = factor(site))) +
    geom_point(size = 3) +  # Scatter plot of the original data
    geom_line(data = new_data_C_pred, aes(y = predicted), color = 'darkgreen', size = 1) +  # Regression line
    labs(title = "Mixed-Effects Regression Line Plot", 
         x = "Silicate (umol/L)", 
         y = "umol/cm^2/hr") +  # Labels
    theme_minimal()

Pred_plot_C_10

### it's working rn with this code but now I only have one site line ? 

Comparing Actual vs Predicted Nutrient Values

  • these plots are meant to compare the actual env data collected from the dilutions to the predicted values based on calculations from the endmember gw data
ActualvsPredicted_Nutrients_edit_Nitrate <- ActualvsPredicted_Nutrients %>% 
  pivot_longer(cols = "Actual":"Predicted", names_to = "Value_Type", values_to = "Value") %>% 
  filter(Parameter=="Nitrate")

ActualvsPredicted_Nutrients_edit_Phosphate <- ActualvsPredicted_Nutrients %>% 
  pivot_longer(cols = "Actual":"Predicted", names_to = "Value_Type", values_to = "Value") %>% 
  filter(Parameter=="Phosphate") %>% 
  filter(Value<1)

#### averaging the values across days per SGD dil to get just one value 
ActualvsPredicted_Nutrients_edit_Nitrate_averaged <- ActualvsPredicted_Nutrients_edit_Nitrate %>% 
  group_by(SGD_number, Site, Parameter, Value_Type) %>% 
  dplyr::summarize(meanN = mean(Value))

ActualvsPredicted_Nutrients_edit_Phosphate_averaged <- ActualvsPredicted_Nutrients_edit_Phosphate %>% 
  group_by(SGD_number, Site, Value_Type) %>% 
  dplyr::summarize(meanN = mean(Value))
  
#########################
## nitrate 
#########################
 ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED <- ActualvsPredicted_Nutrients_edit_Nitrate_averaged %>%
 ggplot(aes(x=(SGD_number+1), 
            y=meanN, 
            color=Value_Type)) + 
  geom_point(size=4) +
 scale_color_manual(values=pnw_palette("Bay", n=2)) +
  facet_wrap(~Site) +
# geom_smooth(method="lm", formula="y~log(x)", color="black") + 
  theme_bw() + 
 coord_trans(x ="log") +
  labs(x = "SGD Percentages",
       y = expression(Nitrate+Nitrite~(mu*mol~L^-1)),
      # title = "Predicted vs Actual Values by Nitrate",
       color = "Nutrient Measurements") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20),
        axis.title.y=element_text(size=20),
        legend.title = element_text(size=18), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED

ggsave(here::here("Outputs", "Results", "Actual_vs_Predicted_Nitrates_byday_AVG.jpg"), 
       width = 10, height=7)

#########################
## phosphate 
#########################
ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED <- ActualvsPredicted_Nutrients_edit_Phosphate_averaged %>% 
 ggplot(aes(x=(SGD_number+1), 
            y=meanN, 
            color=Value_Type)) + 
  geom_point(size=4) +
 scale_color_manual(values=pnw_palette("Bay", n=2)) +
  facet_wrap(~Site) +
# geom_smooth(method="lm", formula="y~log(x)", color="black") + 
  theme_bw() + 
 coord_trans(x ="log") +
  labs(x = "SGD Percentages",
      y = expression(Phosphate~(mu*mol~L^-1)),
      # title = "Predicted vs Actual Values by Phosphate",
       color = "Nutrient Measurements") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20),
        axis.title.y=element_text(size=20),
        legend.title = element_text(size=18), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED

ggsave(here::here("Outputs", "Results", "Actual_vs_Predicted_Phosphate_nooutlier_AVG.jpg"), 
       width = 10, height=7)
  

### patch them together 
p1.1a <- ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.1b <- ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED 

actualVpredicted_NutsPlot <- p1.1a/p1.1b +
  plot_layout(guides = 'collect') + 
  plot_annotation(tag_levels = 'A') #+
  #theme(axis.title = element_text(face = "bold"), 
     #   plot.tag = element_text(size = 12, face = "bold"))

actualVpredicted_NutsPlot

ggsave(here::here("Outputs", "Results", "Final", "ActualVPredicted_NutsPlot.jpg"), 
       width = 12, height=8)


#########################################
## look at day by day measurements 
## because averages are off for first three dilutions for actual data 
#########################################

### NN
NN_byday_C <- ActualvsPredicted_Nutrients_edit_Nitrate %>% 
  filter(Site=="Cabral") %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=Value, 
             color=Value_Type)) + 
  geom_point() +
  facet_wrap(~Day)  + 
  scale_color_manual(values=pnw_palette("Bay", n=2)) +
   coord_trans(x ="log") + 
  labs(title="NN for Cabral")

NN_byday_C

NN_byday_V <- ActualvsPredicted_Nutrients_edit_Nitrate %>% 
  filter(Site=="Varari") %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=Value, 
             color=Value_Type)) + 
  geom_point() +
  facet_wrap(~Day)  + 
  scale_color_manual(values=pnw_palette("Bay", n=2)) +
   coord_trans(x ="log") + 
  labs(title="NN for Varari")

NN_byday_V

### Phos
Phos_byday_C <- ActualvsPredicted_Nutrients_edit_Phosphate %>% 
  filter(Site=="Cabral") %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=Value, 
             color=Value_Type)) + 
  geom_point() +
  facet_wrap(~Day) + 
  scale_color_manual(values=pnw_palette("Bay", n=2)) +
   coord_trans(x ="log") + 
  labs(title="Phosphate for Cabral")

Phos_byday_C

Phos_byday_V <- ActualvsPredicted_Nutrients_edit_Phosphate %>% 
  filter(Site=="Varari") %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=Value, 
             color=Value_Type)) + 
  geom_point() +
  facet_wrap(~Day) + 
  scale_color_manual(values=pnw_palette("Bay", n=2)) +
   coord_trans(x ="log") + 
  labs(title="Phosphate for Varari")

Phos_byday_V

#########################
### Patch 
#########################
NN_byday_C + Phos_byday_C + NN_byday_V + Phos_byday_V

ggsave(here::here("Outputs", "Results", "exploring_actualvspred.jpg"), 
       width = 12, height=7)

check GP and C again colored by site and TA by NEC rates

  • Results show that TA and NEC relationship is opposite of what we expect … so it is likely that nutrients have a stronger effect on the physiological rates than TA
### GP on x and C on y - colored by site 
RespoR_Normalized_Full_withNEC_nutrients2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
  select(!c("phosphate_umolL":"TA_initial")) %>% 
  filter(P_R!="R", P_R!="NP") %>% 
  pivot_wider(names_from = P_R, values_from = umol.cm2.hr)


GP_bothsites_forNEC <- RespoR_Normalized_Full_withNEC_nutrients2 %>% 
  filter(sample_ID!="Varari_Col8_Dil6_Light") %>% 
 # filter(P_R=="C") %>% 
  filter(!is.infinite(GP)) %>% 
  ggplot(aes(x=GP, 
             y=C)) + 
  geom_point(size=4, aes(color=site)) +
  scale_color_manual(values=pnw_palette("Bay", n=2), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "GP",
       y = "C",
       title = "GP by C both sites") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_bothsites_forNEC

## model 
model_GP_C_bothsites <- lmer(GP ~ C*site + (1|new_colonynumber), data = RespoR_Normalized_Full_withNEC_nutrients2 %>% 
                                filter(sample_ID!="Varari_Col8_Dil6_Light") %>% 
                                filter(!is.infinite(GP)))
anova(model_GP_C_bothsites)
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF   DenDF F value Pr(>F)
## C      0.032975 0.032975     1 112.006  1.9160 0.1691
## site   0.000061 0.000061     1  19.805  0.0035 0.9533
## C:site 0.006054 0.006054     1 112.006  0.3517 0.5543
############################
## TA by NEC 
###########################

TA_bothsites_forNEC <- RespoR_Normalized_Full_withNEC_nutrients %>% 
#filter(sample_ID=="Varari_Col8_Dil6_Light") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=TA_initial, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=site)) +
  scale_color_manual(values=pnw_palette("Bay", n=2), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "TA",
       y = "C",
       title = "TA by C both sites") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
TA_bothsites_forNEC

For ECRS:

  • editing the calcification plots specifically for the ECRS talk
###############################
## Calcification 
###############################

C_silicate_ECRS_firebrick4 <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  mutate(site = dplyr::recode(site, 
                      "Varari" = "Site 2",
                      "Cabral" = "Site 1")) %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=5, aes(color=site)) +
  scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
  theme_classic() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  geom_hline(yintercept = 0) +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "C", site == "Varari", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y ~ log(x), color = "firebrick4") +
  geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
                filter(P_R == "C", site == "Cabral", !is.infinite(umol.cm2.hr)), 
              aes(x = silicate_umolL, y = umol.cm2.hr), 
              method = "lm", formula = y ~ log(x), color = "darkgoldenrod2", linetype = "dashed") +
 # geom_smooth(method="lm", formula = "y~log(x)", color="red3", 
           #   method="lm", formula = "y~log(x)", color="deepskyblue4", type="dashed") +
 # geom_smooth(method="lm", formula = "y~log(x)", color="deepskyblue4", type="dashed") +
  labs(x = "Logged Silicate (umolL)",
       y = "Rate (umol CaCO3 g-1 hr-1)",
       title = "Calcification Rates as a Function of Silicate",
       color = "Site") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=20, face="bold"),
        axis.title.y=element_text(size=20, face="bold"),
        legend.title = element_text(size=18, 
                                    face="bold"), 
        legend.text = element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20, 
                                face="bold"))
C_silicate_ECRS_firebrick4

ggsave(here::here("Outputs", "Results", "C_silicate_ECRS.jpg"), 
       width=10, height=7)

EXtra:

  • not really using and not commented but saving as back up
#RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP <- RespoR_Normalized_Full_withNEC_nutrients %>% 
#  filter(site=="Cabral") %>% 
 # filter(P_R=="C") %>% 
#  filter(!is.infinite(umol.cm2.hr))

#min(RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP$umol.cm2.hr)
#max(RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP$umol.cm2.hr)




#RespoR_Normalized_Full_withNEC_nutrients_Varari_GP <- RespoR_Normalized_Full_withNEC_nutrients %>% 
 # filter(site=="Varari") %>% 
 #  filter(P_R=="C") %>% 
 # filter(!is.infinite(umol.cm2.hr))

#min(RespoR_Normalized_Full_withNEC_nutrients_Varari_GP$umol.cm2.hr)
#max(RespoR_Normalized_Full_withNEC_nutrients_Varari_GP$umol.cm2.hr)


### plot 

 #C_silicate_edit <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  #filter(P_R=="C") %>% 
 # filter(!is.infinite(umol.cm2.hr)) %>% 
  #ggplot(aes(x=silicate_umolL, 
          #   y=umol.cm2.hr)) + 
 # geom_point(size=4, aes(color=site)) +
 # scale_color_manual(values=pnw_palette("Bay", n=2)) +
  #theme_classic() + 
  #coord_trans(x ="log") +
  #geom_abline(slope = -0.02219, intercept = 0.15352, color = "black") +
 # geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2)
  #geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>% 
              #  filter(P_R == "C", site == "Varari", !is.infinite(umol.cm2.hr)), 
            #  aes(x = silicate_umolL, y = umol.cm2.hr), 
             # method = "lm", formula = y ~ log(x), color = "red3") 
 # labs(x = "logged Silicate (umolL)",
 #      y = "Rate (umol CaCO3 g-1 hr-1)",
     #  title = "Calcification Rates as a Function of Silicate",
    #   color = "Site") +
  # theme(axis.text.x=element_text(size=18), 
      #  axis.text.y=element_text(size=18), 
     #   axis.title.x=element_text(size=20),
     #   axis.title.y=element_text(size=20),
     ##   legend.title = element_text(size=18), 
      #  legend.text = element_text(size=18), 
      #  plot.title=element_text(hjust=0.5, ## centers title 
                  #              size=20))
#C_silicate_edit





##############################################################
##############################################################

###############################
### trying to use marginal and conditional values to plot 
###############################

### mixed model 
#C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
                                         #   filter(P_R=="C") %>% 
                                          #  filter(!is.infinite(umol.cm2.hr)))

#C_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
 # filter(P_R=="C") %>% 
 # filter(!is.infinite(umol.cm2.hr)) %>% 
 # mutate(fit.m = predict(C_bothsites_silicate_linearmodel),
        # fit.c = predict(C_bothsites_silicate_linearmodel)) #+ 
 # mutate(resid = resid(C_bothsites_silicate_linearmodel)) ### getting an error here


#C_silicate2.1 <- C_silicate2 %>% 
 # ggplot(aes(x=silicate_umolL, 
 #            y=umol.cm2.hr)) + 
#  geom_point(pch = 16, col = "grey") +
 #geom_line(aes(y = fit.c), col = 1, size = 2) 
  
##C_silicate2.1  
  
 # geom_point(size=4, aes(color=site)) +
 # scale_color_manual(values=pnw_palette("Bay", n=2)) +
  #theme_bw() + 
 # coord_trans(x ="log") +
 # theme(strip.background = element_rect(fill = "white")) +
#  geom_hline(yintercept = 0) +
 # geom_line(aes(y = fit.c), col = 1, size = 2) +
#  coord_cartesian(ylim = c(-40, 100)) +
  #geom_smooth(method="lm", formula = "y~log(x)", color="black") +
#  labs(x = "logged Silicate (umolL)",
  #     y = "Rate (umol O2 cm2 hr-1)",
   #    title = "Calcification Rates as a Function of Silicate",
   #    color = "Site") +
 #  theme(axis.text.x=element_text(size=18), 
     #   axis.text.y=element_text(size=18), 
     #   axis.title.x=element_text(size=18),
    #    axis.title.y=element_text(size=18), 
    #    plot.title=element_text(hjust=0.5, ## centers title 
                          #      size=20))
#C_silicate2.1#

#ggsave(here::here("Outputs", "Results", "C_silicate.jpg"))




###############################
## for Calcification 
###############################
#RespoR_Normalized_Full_withNEC_nutrients2<- RespoR_Normalized_Full_withNEC_nutrients %>%
  #filter(P_R=="C") %>% 
#  filter(!is.infinite(umol.cm2.hr)) %>% 
#  arrange(RespoR_Normalized_Full_withNEC_nutrients2$umol.cm2.hr) 

#C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), #data=RespoR_Normalized_Full_withNEC_nutrients2)

### ^^ look for slope and intercept 
#summary(C_bothsites_silicate_linearmodel)


### use predict function to get info for confidence intervals - use cbind to create a new df and use that 

#predC <- predict(C_bothsites_silicate_linearmodel) # Gets the predicted values from the regression lines in the ANOVA

#graphdata_C <- cbind(RespoR_Normalized_Full_withNEC_nutrients2, predC) # attaches those predictions to the dataset

#graphdata_C <- graphdata_C %>% 
  #arrange(graphdata_C$silicate_umolL)


# We'll plot the original data as points, and the regression predictions as lines 
#C_silicate_pred <- new_data %>% 
 # ggplot(aes(x=silicate_umolL, 
    #         y=umol.cm2.hr, 
       #      color=site)) +
 # geom_point(size=4) + 
 # geom_line(aes(y=(predC))) +
#  geom_line(data = new_data, aes(x = log(silicate_umolL), y = predictions, color = site)) +
#  scale_color_manual(values=pnw_palette("Bay", n=2)) +
#  theme_classic() + 
#  coord_trans(x ="log") +
#  labs(x = "logged Silicate (umolL)",
    #   y = "Rate (umol CaCO3 g-1 hr-1)",
    #   title = "Calcification Rates as a Function of Silicate",
    #   color = "Site") +
 #  theme(axis.text.x=element_text(size=18), 
    #    axis.text.y=element_text(size=18), 
    #    axis.title.x=element_text(size=20),
    #    axis.title.y=element_text(size=20),
    #    legend.title = element_text(size=18), 
    #    legend.text = element_text(size=18), 
     #   plot.title=element_text(hjust=0.5, ## centers title 
                           #     size=20))
#C_silicate_pred


###############################
## for Gross Photosynthesis  
###############################
#RespoR_Normalized_Full_withNEC_nutrients3 <- RespoR_Normalized_Full_withNEC_nutrients %>% 
 # # filter(P_R=="GP") %>% 
 # filter(!is.infinite(umol.cm2.hr))

#GP_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), #data=RespoR_Normalized_Full_withNEC_nutrients3)

### ^^ look for slope and intercept 
#summary(GP_bothsites_silicate_linearmodel)

############################
### attemot 1 -- something is off with the pred function, this is for a simpler model, mine is too complex 
############################

#Cpred <- predict(C_bothsites_silicate_linearmodel)

#graph_Cpred_data <- cbind(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted, Cpred) 

#Cpred_plot <- ggplot(data=graph_Cpred_data,
    #   aes(x=log(silicate_umolL), 
     #      y=umol.cm2.hr, 
     #      color=site)) +
 # geom_point(size=4) + 
 # geom_line(aes(y=Cpred)) +
 # labs(x="silicate_umolL", y="umol.cm2.hr", fill="site") +
 # theme_bw() 

#Cpred_plot

#Cpred_plot2 <- ggplot(data=graph_Cpred_data,
#       aes(x=log(silicate_umolL), 
 #          y=Cpred, 
 #          color=site)) +
#  geom_point(size=4) + 
#  geom_line(aes(y=Cpred)) +
 # labs(x="silicate_umolL", y="Cpred", fill="site") +
 # theme_bw() 


#Cpred_plot2






###########################
## get basic correlations with p-values 
##########################
#res.cor.Varari <- correlate(corr_matrix_Varari, method = "spearman", diagonal = 1)

###########################
## edit correlation matrix and plot 
##########################
#res.cor.Varari %>% 
  #rearrange(method = "MDS", absolute = FALSE) %>% 
 # shave(upper=TRUE) %>% 
 # filter(!Salinity = 1) %>% 
 # rplot() ## make correlologram 
  

###########################
#res.cor.Varari %>% 
  #filter(!"Salinity" = 1)


 # rearrange(method = "MDS", absolute = FALSE) %>% 
 # shave(upper=TRUE) %>% 
  
##########################

Make correlation matrix for nutrients

## data 
corr_matrix_envparams <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="R") %>% 
  mutate(Salinity=salinity) %>% 
  select(!c("salinity", "umol.cm2.hr", "SGD_number", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) ## removed ammonia for rn because have less than values and causing NAs

##visualize it 
ggcorr(corr_matrix_envparams, method = c("everything", "pearson")) 

###############################
### Dani B's methods #######
cor_mat <- round(cor(corr_matrix_envparams),4)
head(cor_mat)
##                 phosphate_umolL silicate_umolL NN_umolL  temp_C  new_pH
## phosphate_umolL          1.0000         0.2483   0.2789  0.1626 -0.0260
## silicate_umolL           0.2483         1.0000   0.1296 -0.1146 -0.1927
## NN_umolL                 0.2789         0.1296   1.0000  0.1277 -0.0638
## temp_C                   0.1626        -0.1146   0.1277  1.0000  0.2110
## new_pH                  -0.0260        -0.1927  -0.0638  0.2110  1.0000
## TA_initial                   NA             NA       NA      NA      NA
##                 TA_initial log_silicate log_silicate2 Salinity
## phosphate_umolL         NA       0.2378        0.2483  -0.2635
## silicate_umolL          NA       0.9590        0.9887  -0.6914
## NN_umolL                NA       0.1640        0.1465   0.2955
## temp_C                  NA      -0.1790       -0.1415  -0.1727
## new_pH                  NA      -0.2150       -0.2027  -0.0230
## TA_initial               1           NA            NA       NA
#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
melted_cormat <- melt(cor_mat)
head(melted_cormat)
##              Var1            Var2   value
## 1 phosphate_umolL phosphate_umolL  1.0000
## 2  silicate_umolL phosphate_umolL  0.2483
## 3        NN_umolL phosphate_umolL  0.2789
## 4          temp_C phosphate_umolL  0.1626
## 5          new_pH phosphate_umolL -0.0260
## 6      TA_initial phosphate_umolL      NA
#visulaize the correlation matrix in general
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile()

# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
  cormat[upper.tri(cor_mat)] <- NA
  return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
  cor_mat[lower.tri(cor_mat)]<- NA
  return(cor_mat)
}

#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
upper_tri
##                 phosphate_umolL silicate_umolL NN_umolL  temp_C  new_pH
## phosphate_umolL               1         0.2483   0.2789  0.1626 -0.0260
## silicate_umolL               NA         1.0000   0.1296 -0.1146 -0.1927
## NN_umolL                     NA             NA   1.0000  0.1277 -0.0638
## temp_C                       NA             NA       NA  1.0000  0.2110
## new_pH                       NA             NA       NA      NA  1.0000
## TA_initial                   NA             NA       NA      NA      NA
## log_silicate                 NA             NA       NA      NA      NA
## log_silicate2                NA             NA       NA      NA      NA
## Salinity                     NA             NA       NA      NA      NA
##                 TA_initial log_silicate log_silicate2 Salinity
## phosphate_umolL         NA       0.2378        0.2483  -0.2635
## silicate_umolL          NA       0.9590        0.9887  -0.6914
## NN_umolL                NA       0.1640        0.1465   0.2955
## temp_C                  NA      -0.1790       -0.1415  -0.1727
## new_pH                  NA      -0.2150       -0.2027  -0.0230
## TA_initial               1           NA            NA       NA
## log_silicate            NA       1.0000        0.9899  -0.5874
## log_silicate2           NA           NA        1.0000  -0.6489
## Salinity                NA           NA            NA   1.0000
#melt the correlation matrix
#melt the correlation data and drop the rows with NA values 
melted_cormat <- melt(upper_tri, na.rm = TRUE)

#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()

melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))

# Create a ggheatmap with basic characteristics, etc. 
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  theme_minimal()+ # minimal theme
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()


# Print the heatmap
print(ggheatmap)

#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color 
corr_matrix_SGD_params <- ggheatmap + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 6) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 18, face="bold", color="black"),
        axis.text.y = element_text(size = 18, face="bold", color="black"),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(0.8, 0),
        legend.title = element_text(size = 18, face="bold", color="black"),
        legend.text = element_text(size = 20, face="bold", color="black"),
        legend.position = c(0.48, 0.75),
        legend.direction = "horizontal") +
  guides(fill = guide_colorbar(barwidth = 12, barheight = 2, 
                               title.position = "top", title.hjust = 0.5, title.vjust = 1.0))

corr_matrix_SGD_params

Figure 2:

############################################
### correlation matrix for Varari 
############################################

## data 
corr_matrix_envparams_Varari <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  mutate(N_P = NN_umolL/phosphate_umolL) %>% 
  filter(P_R=="R") %>% 
  mutate(Salinity=salinity, 
         TA = TA_initial, 
         pH = new_pH, 
         Nitrate_Nitrite_umolL = NN_umolL, 
         Phosphate_umolL = phosphate_umolL,
         Silicate_umolL = silicate_umolL, 
         Ammonia_umolL = ammonia_umolL) %>% 
  filter(site=="Varari") %>% 
  select(!c("salinity", "silicate_umolL", "TA_initial", "Ammonia_umolL", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) 

##visualize it 
ggcorr(corr_matrix_envparams_Varari, method = c("everything", "pearson")) 

###############################
cor_mat <- round(cor(corr_matrix_envparams_Varari),4)
#head(cor_mat)

###############################
# Set the diagonal elements to NA
#diag(cor_mat) <- NA
###############################

###############################
# Melt the correlation matrix
melted_cormat <- melt(cor_mat, na.rm = TRUE)
###############################

#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
#melted_cormat <- melt(cor_mat)
#head(melted_cormat)

#visualize the correlation matrix in general
#ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
#  geom_tile()

# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
  cormat[upper.tri(cor_mat)] <- NA
  return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
  cor_mat[lower.tri(cor_mat)]<- NA
  return(cor_mat)
}

#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
#upper_tri

#melt the correlation matrix
#melt the correlation data and drop the rows with NA values 
melted_cormat <- melt(upper_tri, na.rm = TRUE)

#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson Correlation") +
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()

melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))

# Create a ggheatmap with basic characteristics, etc. 
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson Correlation") +
  theme_minimal()+ # minimal theme
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 10, hjust = 1))+
  coord_fixed() 

diag(melted_cormat) <- NA

# Print the heatmap
#print(ggheatmap)

#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color 
corr_matrix_SGD_params_Varari <- ggheatmap + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 12, face="bold", color="black"),
        axis.text.y = element_text(size = 12, face="bold", color="black"),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(0.8, 0),
        legend.title = element_text(size = 10, color="black"),
        legend.text = element_text(size = 10, color="black"),
        legend.position = c(0.48, 0.75),
        legend.direction = "horizontal", 
        title= element_text(size = 12, hjust=0.5, face="bold", color="black")) +
  labs(title="Pearson Correlation Matrix for Environmental Parameters at Varari") + 
  guides(fill = guide_colorbar(barwidth = 8, barheight = 1, 
                               title.position = "top", title.hjust = 0.5, title.vjust = 1.0))

corr_matrix_SGD_params_Varari

ggsave(here::here("Outputs", "Results", "Varari_correlation_matrix.jpg"), 
       width=10, height=7)

Figure 2a:

############################################
### correlation matrix for Cabral 
############################################

## data 
corr_matrix_envparams_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  mutate(N_P = NN_umolL/phosphate_umolL) %>% 
  filter(P_R=="R") %>% 
  mutate(Salinity=salinity, 
         TA = TA_initial, 
         pH = new_pH, 
         Nitrate_Nitrite_umolL = NN_umolL, 
         Phosphate_umolL = phosphate_umolL,
         Silicate_umolL = silicate_umolL, 
         Ammonia_umolL = ammonia_umolL) %>% 
  filter(site=="Cabral") %>% 
  select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL", "Ammonia_umolL")) 

##visualize it 
ggcorr(corr_matrix_envparams_Cabral, method = c("everything", "pearson")) 

###############################
cor_mat <- round(cor(corr_matrix_envparams_Cabral),4)
#head(cor_mat)

#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
melted_cormat <- melt(cor_mat)
#head(melted_cormat)


#visualize the correlation matrix in general
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile()

# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
  cormat[upper.tri(cor_mat)] <- NA
  return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
  cor_mat[lower.tri(cor_mat)]<- NA
  return(cor_mat)
}

#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
#upper_tri

#melt the correlation matrix
#melt the correlation data and drop the rows with NA values 
melted_cormat <- melt(upper_tri, na.rm = TRUE)

#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson Correlation") +
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()

melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))

# Create a ggheatmap with basic characteristics, etc. 
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson Correlation") +
  theme_minimal()+ # minimal theme
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 10, hjust = 1))+
  coord_fixed()

# Print the heatmap
#print(ggheatmap)

#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color 
corr_matrix_SGD_params_Cabral <- ggheatmap + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 12, face="bold", color="black"),
        axis.text.y = element_text(size = 12, face="bold", color="black"),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(0.8, 0),
        legend.title = element_text(size = 10, color="black"),
        legend.text = element_text(size = 10, color="black"),
        legend.position = c(0.48, 0.75),
        legend.direction = "horizontal", 
        title= element_text(size = 12, hjust=0.5, face="bold", color="black")) +
  labs(title="Pearson Correlation Matrix for Environmental Parameters at Cabral") + 
  guides(fill = guide_colorbar(barwidth = 8, barheight = 1, 
                               title.position = "top", title.hjust = 0.5, title.vjust = 1.0))

corr_matrix_SGD_params_Cabral

ggsave(here::here("Outputs", "Results", "Cabral_correlation_matrix.jpg"), 
       width=10, height=7)

Results of back calculating

Essentially found that the predicted values were too low to be useful for

take average of Varari Days 2-4 to fill in for ambient day 1

  • Did not end up doing, May 15th, 2025
C_only <- RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(!is.infinite(umol.cm2.hr)) %>%
  filter(P_R=="C") %>% 
  mutate(new_colonynumber=as.numeric(new_colonynumber))


-18.427160
## [1] -18.42716
-57.911680
## [1] -57.91168
-63.030418
## [1] -63.03042
(-18.427160+-57.911680+-63.030418)/3
## [1] -46.45642
## REPLACE 146 FROM AMBIETN DAY 1 WITH -46.45642 WHICH IS AVERAGE OF ALL 3 OTHER DAYS FOR DELTA TA AT SGD NUMBER 1, AMB SW